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We report on an all-sky search with the LIGO detectors for periodic gravitational waves in the 
frequency range 50- 1000 Hz and with the frequency's time derivative in the range —1 x 10~* Hz s~^ 
to zero. Data from the fourth LIGO science run (S4) have been used in this search. Three different 
semi-coherent methods of transforming and summing strain power from Short Fourier Transforms 
(SFTs) of the calibrated data have been used. The first, known as "StackSlide" , averages normalized 
power from each SFT. A "weighted Hough" scheme is also developed and used, and which also 
allows for a multi-interferometer search. The third method, known as "PowerFlux", is a variant 
of the StackSlide method in which the power is weighted before summing. In both the weighted 
Hough and PowerFlux methods, the weights are chosen according to the noise and detector antenna- 
pattern to maximize the signal-to-noise ratio. The respective advantages and disadvantages of these 
methods are discussed. Observing no evidence of periodic gravitational radiation, we report upper 
limits; we interpret these as limits on this radiation from isolated rotating neutron stars. The best 
population-based upper limit with 95% confidence on the gravitational-wave strain amplitude, found 
for simulated sources distributed isotropically across the sky and with isotropically distributed spin- 
axes, IS 4.28 X 10"^" (near 140 Hz). Strict upper limits are also obtained for small patches on the 
sky for best-case and worst-case inclinations of the spin axes. 

PACS numbers: 04.80.Nn, 95.55.Ym, 97.60.Gb, 07.05.Kf 



I. INTRODUCTION 

We report on a search with the LIGO (Laser Interfer- 
ometer Gravitational-wave Observatory) detectors [TJ [5] 
for periodic gravitational waves in the frequency range 
50 - 1000 Hz and with the frequency's time derivative in 
the range —1 x 10^^ Hz s^^ to zero. The search is carried 
out over the entire sky using data from the fourth LIGO 
science run (S4). Isolated rotating neutron stars in our 
galaxy are the prime target. 

Using data from earlier science runs, the LIGO Sci- 
entific Collaboration (LSC) has previously reported on 
searches for periodic gravitational radiation, using a long- 
period coherent method to target known pulsars [3 HI E] , 
using a short-period coherent method to target Scorpius 
X-1 in selected bands and search the entire sky in the 
160.0-728.8 Hz band [B], and using a long-period semi- 
coherent method to search the entire sky in the 200- 
400 Hz band [7]. Einstein@Homc, a distributed home 
computing effort running under the BOINC architecture 
has also been searching the entire sky using a coher- 
ent first stage, followed by a simple coincidence stage [S]. 
In comparison, this paper: 1) examines more sensitive 
data; 2) searches over a larger range in frequency and its 
derivative; and 3) uses three alternative semi-coherent 
methods for summing measured strain powers to detect 
excess power from a continuous gravitational- wave signal. 

The first purpose of this paper is to present results from 
our search for periodic gravitational waves in the S4 data. 
Over the LIGO frequency band of sensitivity, the S4 all- 
sky upper limits presented here are approximately an or- 
der of magnitude better than published previously from 
earlier science runs |B] |7] . After following up on outliers 
in the data, we find that no candidates survive, and thus 
report upper limits. These are interpreted as limits on 
radiation from rotating neutron stars, which can be ex- 
pressed as functions of the star's ellipticity and distance, 
allowing for an astrophysical interpretation. The best 
population-based upper limit with 95% confidence on the 
gravitational- wave strain amplitude, found for simulated 



sources distributed isotropically across the sky and with 
isotropically distributed spin-axes, is 4.28 x 10^^^ (near 
140 Hz). Strict upper limits are also obtained for small 
patches on the sky for best-case and worst-case inclina- 
tions of the spin axes. 

The second purpose of this paper, along with the pre- 
vious coherent [6 and semi-coherent 7^ papers, is to lay 
the foundation for the methods that will be used in fu- 
ture searches. It is well known that the search for periodic 
gravitational waves is computationally bound; to obtain 
optimal results will require a hierarchical approach that 
uses coherent and semi-coherent stages (TUl [HI EH US] . A 
fifth science run (S5), which started in November 2005, is 
generating data at initial LIGO's design sensitivity. We 
plan to search this data using the best methods possible, 
based on what is learned from this and previous analyses. 

In the three methods considered here, one searches 
for cumulative excess power from a hypothetical periodic 
gravitational wave signal by examining successive spec- 
tral estimates based on Short Fourier Transforms (SFTs) 
of the calibrated detector strain data channel, taking into 
account the Doppler modulations of detected frequency 
due to the Earth's rotational and orbital motion with 
respect to the Solar System Barycenter (SSB), and the 
time derivative of the frequency intrinsic to the source. 
The simplest method presented, known as "StackSlide" 
[121 [131 El [IS], averages normalized power from each 
SFT. In the Hough method reported previously [71 [TO], 
referred to here as "standard Hough" , the sum is of bi- 
nary zeroes or ones, where an SFT contributes unity if 
the power exceeds a normalized power threshold. In this 
paper a "weighted Hough" scheme, henceforth also re- 
ferred to as "Hough", has been developed and is simi- 
lar to that described in Ref. [TB]. This scheme also al- 
lows for a multi-interferometer search. The third method, 
known as "PowerFlux" [T7] , is a variant of the StackSlide 
method in which the power is weighted before summing. 
In both the weighted Hough and PowerFlux methods, 
the weights are chosen according to the noise and de- 
tector antenna pattern to maximize the signal-to-noise 
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ratio. 

The Hough method is computationally faster and more 
robust against large transient power artifacts, but is 
slightly less sensitive than StackSlide for stationary data 
[TJ [TS]. The PowerFlux method is found in most fre- 
quency ranges to have better detection efficiency than 
the StackSlide and Hough methods, the exceptions oc- 
curring in bands with large non-stationary artifacts, for 
which the Hough method proves more robust. However, 
the StackSlide and Hough methods can be made more 
sensitive by starting with the maximum likelihood statis- 
tic (known as the J^-statistic |S1[TDJ[TH]) rather than SFT 
power as the input data, though this improvement comes 
with increased computational cost. The trade-offs among 
the methods means that each could play a role in our fu- 
ture searches. 

In brief, this paper makes several important contribu- 
tions. It sets the best all-sky upper limits on periodic 
gravitational waves to date, and shows that these limits 
are becoming astrophysically interesting. It also intro- 
duces methods that are crucial to the development of 
our future searches. 

This paper is organized as follows: Section |II] briefly 
describes the LIGO interferometers, focusing on improve- 
ments made for the S4 data run, and discusses the sen- 
sitivity and relevant detector artifacts. Section |ITT] pre- 
cisely defines the waveforms we seek and the associated 



assumptions we have made. Section IV gives a detailed 
description of the three analysis methods used and sum- 
marizes their similarities and differences, while Section 
[V| gives the details of their implementations and the 



pipelines used. Section VI discusses the validation of the 
software and, as an end-to-end test, shows the detection 
of simulated pulsar signals injected into the data stream 
at the hardware level. Section IVIII describes the search 



results, and Section VIII compares the results from the 
three respective methods. Section |IX| concludes with a 
summary of the results, their astrophysical implications, 
and future plans. 



II. THE LIGO DETECTOR NETWORK AND 
THE S4 SCIENCE RUN 

The LIGO detector network consists of a 4-km inter- 
ferometer in Livingston Louisiana (called LI) and two 
interferometers in Hanford Washington, one 4-km and 
another 2-km (HI and H2, respectively). 

The data analyzed in this paper were produced dur- 
ing LIGO's 29.5-day fourth science run (S4) [19]. This 
run started at noon Central Standard Time (GST) on 
February 22 and ended at midnight GST on March 23, 
2005. During the run, all three LIGO detectors had dis- 
placement spectral amplitudes near 2.5 x 10^^^ m Hz^^^^ 
in their most sensitive frequency band near 150 Hz. In 
units of gravitational-wave strain amplitude, the sensi- 
tivity of H2 is roughly a factor of two worse than that 
of HI and LI over much of the search band. The typical 



strain sensitivities in this run were within a factor of two 
of the design goals. Figure [l] shows representative strain 
spectral noise densities for the three interferometers dur- 
ing the run. As discussed in Section |V] below, however, 
non-stationarity of the noise was significant. 

Ghanges to the interferometers before the S4 run in- 
cluded the following improvements |19j : 

• Installation of active seismic isolation of support 
structures at Livingston to cope with high anthro- 
pogenic ground motion in the 1-3 Hz band. 

• Thermal compensation with a GO2 laser of mirrors 
subject to thermal lensing from the primary laser 
beam to a greater or lesser degree than expected. 

• Replacement of a synthesized radio frequency oscil- 
lator for phase modulation with a crystal oscillator 
before S4 began (HI) and mid- way through the S4 
run (LI), reducing noise substantially above 1000 
Hz and eliminating a comb of ~ 37 Hz lines. (The 
crystal oscillator replacement for H2 occurred after 
the S4 run.) 

• Lower-noise mirror-actuation electronics (HI, H2, 
& LI). 

• Higher-bandwidth laser frequency stabilization 
(HI, H2, & LI) and intensity stabilization (HI & 
LI). 

• Installation of radiation pressure actuation of mir- 
rors for calibration validation (HI). 

• Gommissioning of complete alignment control sys- 
tem for the LI interferometer (already implemented 
for HI & H2 in S3 run). 

• Refurbishment of lasers and installation of photo- 
diodes and electronics to permit interferometer op- 
eration with increased laser power (HI, H2, & LI). 

• Mitigation of electromagnetic interference (HI, H2, 
& LI) and acoustic interference (LI). 

The data were acquired and digitized at a rate of 16384 
Hz. Data acquisition was periodically interrupted by dis- 
turbances such as seismic transients, reducing the net 
running time of the interferometers. The resulting duty 
factors for the interferometers were 81% for HI and H2, 
and 74% for LI. While the HI and H2 duty factors were 
somewhat higher than those in previous science runs, the 
LI duty factor was dramatically higher than the ~40% 
typical of the past, thanks to the increased stability from 
the installation of the active seismic isolation system at 
Livingston. 



III. SIGNAL WAVEFORMS 

The general form of a gravitational-wave signal is de- 
scribed in terms of two orthogonal transverse polariza- 
tions defined as with waveform h^{t) and "x" with 
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FIG. 1: Median amplitude strain noise spectral densities from 
the three LIGO interferometers during the S4 run, along with 
the Initial LIGO design sensitivity goal. 



waveform hx(t). The calibrated response seen by an in- 
terferometric gravitational-wave detector is then [18] 

h{t) = F+{t, a, S, ^)h+{t) + Fx {t, a, 6, i^)hx (t), (1) 

where t is time in the detector frame, a is the source 
right ascension, 6 is the source declination, ip is the po- 
larization angle of the wave, and -F+,x are the detector 
antenna pattern functions for the two orthogonal polar- 
izations. For periodic (nearly pure sinusoidal) gravita- 
tional waves, which in general are elliptically polarized, 
the individual components /i+,x have the form 



h+{t) 
hx{t) 



A+ cos<i>(t), 
ylxsin$(t), 



(2) 
(3) 



where Aj^ and are the amplitudes of the two polariza- 
tions, and <^{t) is the phase of the signal at the detector. 
(One can also define the initial phase of the signal, <i>o, 
but in this paper it can be taken to be an unknown and 
irrelevant constant). 

For an isolated quadrupolar gravitational-wave emit- 
ter, characterized by a rotating triaxial ellipsoid mass 
distribution, the amplitudes and A^ are related to 
the inclination angle of the source, t, and the wave am- 
plitude, /iQ, by: 



A, 



1 



cos b 



Ax = /iQCOst, 



(4) 
(5) 



where i is the angle of its spin axis with respect to the 
line of sight between source and detector. For such a star, 
the gravitational-wave frequency, /, is twice the rotation 
frequency, v, and the amplitude is given by 



(6) 



Here d is the distance to the star, / is the principal mo- 
ment of inertia with respect to its spin axis, and e is the 
equatorial ellipticity of the star [TB] . Assuming that all of 
the frequency's derivative, /, is due to emission of grav- 
itational radiation and that / takes the canonical value 
10^8 kgm^, we can relate e to / and / and use Eq. ^ to 
obtain 



h^^ = 4.54 X 10 
by eliminating e, or 



,24^ IkpcN^ ( 250 yr V 

-fim) ' 



esd = 7.63 X 10 



-5 



10-1" Hz s-i 



100 Hz\ ^ 



, (8) 



by eliminating d. These are referred to, respectively, 
as the spin-down limits on strain and ellipticity. (See 
Eqs. (8), (9), and (19) of ^ for more details of the deriva- 
tion.) 

Note that the methods used in this paper are sen- 
sitive to periodic signals from any type of isolated 
gravitational- wave source (e.g., freely precessing or os- 
cillating neutron stars as well as triaxial ones), though 
we present upper limits in terms of and e. Be- 
cause we use semi-coherent methods, only the instan- 
taneous signal frequency in the detector reference frame, 
2n f{t) ~ d<^{t)/dt, needs to be calculated. In the de- 
tector reference frame this can, to a very good approx- 
imation, be related to the instantaneous SSB-frame fre- 
quency f{t) by ^ 



fit) - fit) - m 



v(t) • n 



(9) 



where v{t) is the detector's velocity with respect to the 
SSB frame, and n is the unit- vector corresponding to the 
sky-location of the source. In this analysis, we search 
for f(t) signals well described by a nominal frequency fo 
at the start of the S4 run tQ and a constant first time 
derivative /, such that 



.fit) = fo + fit-to) 



(10) 



These equations ignore corrections to the time interval 
t — to at the detector compared with that at the SSB and 
relativistic corrections. These corrections are negligible 
for the one month semi-coherent searches described here, 
though the LSC Algorithm Library (LAL) code [2U] used 
by our searches does provide routines that make all the 
corrections needed to provide a timing accuracy of 3 fj,s. 
(The LAL code also can calculate f{t) for signals arriving 
from periodic sources in binary systems. Including un- 
known orbital parameters in the search, however, would 
greatly increase the computational cost or require new 
methods beyond the scope of this article.) 
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FIG. 2: An illustration of the discrete frequency bins of the 
Short Fourier Transform (SFTs) of the data are shown verti- 
cally, with the discrete start times of the SFTs shown hori- 
zontally. The dark pixels represent a signal in the data. Its 
frequency changes with time due to Doppler shifts and in- 
trinsic evolution of the source. By sliding the frequency bins, 
the power from a source can be lined up and summed after 
appropriate weighting or transformation. This is, in essence, 
the starting point for all of the semi-coherent search meth- 
ods presented here, though the actual implementations differ 
significantly. 



IV. OVERVIEW OF THE METHODS 

A. Similarities and Differences 

The three different analysis methods presented here 
have many features in common, but also have important 
differences, both major and minor. In this Section we 
give a brief overview of the methods. 



1. The parameter space 

All three methods are based on summing measures of 
strain power from many SFTs that have been created 
from 30-minute intervals of calibrated strain data. Each 
method also corrects explicitly for sky-position depen- 
dent Doppler modulations of the apparent source fre- 
quency due to the Earth's rotation and its orbital motion 
around the SSB, and the frequency's time derivative, in- 
trinsic to the source (see Fig.|2]). This requires a search 
in a four-dimensional parameter space; a template in the 
space refers to a set of values: A = {/o, /, a, (5}. The 
third method, PowerFlux, also searches explicitly over 
polarization angle, so that A — {/o, /, a, (5, 

All three methods search for initial frequency /o in the 
range 50-1000 Hz with a uniform grid spacing equal to 
the size of an SFT frequency bin, 

5/ = = 5.556 X 10"* Hz . (11) 

^ coh 

where Tcoh is the time-baseline of each SFT. The range 
of /o is determined by the noise curves of the interferom- 
eters, likely detectable source frequencies [21], and limi- 
tations due to the increasing computational cost at high 
frequencies. 



The range of / values searched is [—1 x 10 ^, 0] Hz s ^ 
for the StackSlide and PowerFlux methods and [—2.2 x 
10"^, 0] Hz s"^ for the Hough method. The ranges of / 
are determined by the computational cost, as well as by 
the low probability of finding an object with |/| higher 
than the values searched — in other words, the ranges of / 
are narrow enough to complete the search in a reasonable 
amount of time, yet wide enough to include likely sig- 
nals. All known isolated pulsars spin down more slowly 
than the two values of |/|„ax used here, and as seen in 
the results section, the ellipticity required for higher |/| 
is improbably high for a source losing rotational energy 
primarily via gravitational radiation at low frequencies. 
A small number of isolated pulsars in globular clusters 
exhibit slight spin-up, believed to arise from acceleration 
in the Earth's direction; such spin-up values have magni- 
tudes small enough to be detectable with the zero-spin- 
down templates used in these searches, given a strong 
enough signal. The parameter ranges correspond to a 
minimum spin-down timescale //|4/| (the gravitational- 
wave spin-down age) of 40 years for a source emitting at 
50 Hz and 800 years for a source at 1000 Hz. Since for 
known pulsars [22] this characteristic timescale is at least 
hundreds of years for frequencies on the low end of our 
range and tens of millions of years for frequencies on the 
high end, we see again that the ranges of |/| are wide 
enough to include sources from this population. 

As discussed in our previous reports [SI |7] , the number 
of sky points that must be searched grows quadratically 
with the frequency /o, ranging here from about five thou- 
sand at 50 Hz to about two million at 1000 Hz. All three 
methods use nearly isotropic grids which cover the en- 
tire sky. The PowerFlux search also divides the sky into 
regions according to susceptibility to stationary instru- 
mental line artifacts. Sky grid and spin-down spacings 
and other details are provided below. 

2. Upper limits 

While the parameter space searched is similar for the 
three methods, there are important differences in the way 
upper limits are set. StackSlide and Hough both set 
population-based frequentist limits on by carrying out 
Monte Carlo simulations of a random population of pul- 
sar sources distributed uniformly over the sky and with 
isotropically distributed spin-axes. PowerFlux sets strict 
frequentist limits on circular and linear polarization am- 
plitudes /iCi'-'=-i™it .^^^ /jLin-iimit^ ^j^j^j^ correspond to 

limits on most and least favorable pulsar inclinations, 
respectively. The limits are placed separately on tiny 
patches of the sky, with the highest strain upper limits 
presented here. In this context "strict" means that, re- 
gardless of its polarization angle ■0 or inclination angle t, 
regardless of its sky location (within fiducial regions dis- 
cussed below) , and regardless of its frequency value and 
spin-down within the frequency and spin-down step sizes 
of the search template, an isolated pulsar of true strain 
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amplitude ho = 2/ig 



-limit 



would have yielded a higher 



measured amplitude than what we measure, in at least 
95% of independent observations. The circular polariza- 
tion limits /jCirc-iimit j^ppiy Q^iy iiiQ most favorable 

inclinations (t ~ 0, tt), regardless of sky location and 
regardless of frequency and spin-down, as above. 

Due to these different upper limit setting methods, 
sharp instrumental lines are also handled differently. 
StackSlide and Hough carry out removal of known instru- 
mental lines of varying widths in individual SFTs. The 
measured powers in those bins are replaced with random 
noise generated to mimic the noise observed in neigh- 
boring bins. This line cleaning technique can lead to a 
true signal being missed because its apparent frequency 
may coincide with an instrumental line for a large num- 
ber of SFTs. However, population- averaged upper limits 
are determined self-consistently to include loss of detec- 
tion efficiency due to line removal, by using Monte Carlo 
simulations. 

Since its limits are intended to be strict, that is, valid 
for any source inclination and for any source location 
within its fiducial area, PowerFlux must handle instru- 
mental lines differently. Single-bin lines are flagged dur- 
ing data preparation so that when searching for a partic- 
ular source an individual SFT bin power is ignored when 
it coincides with the source's apparent frequency. If more 
than 80% of otherwise eligible bins are excluded for this 
reason, no attempt is made to set a limit on strain power 
from that source. In practice, however, the 80% cutoff is 
not used because we have found that all such sources lie 
in certain unfavorable regions of the sky, which we call 
"skybands" and which we exclude when setting upper 
limits. These skybands depend on source frequency and 
its derivative, as described in Sec. |VD4 



3. Data Preparation 

Other differences among the methods concern the data 
windowing and filtering used in computing Fourier trans- 
forms and concern the noise estimation. StackSlide and 
Hough apply high pass filters to the data above 40Hz, in 
addition to the filter used to produce the calibrated data 
stream, and use Tukey windowing. PowerFlux applies no 
additional filtering and uses Hann windowing with 50% 
overlap between adjacent SFT's. StackSlide and Hough 
use median-based noise fioor tracking [23l EH El]. In 
contrast, Powerfiux uses a time-frequency decomposition. 
Both of these noise estimation methods are described in 
Sec. El 

The raw, uncalibrated data channels containing the 
strain measurements from the three interferometers are 
converted to a calibrated "h{t)" data stream, following 
the procedure described in [55], using calibration refer- 
ence functions described in [27]. SFTs are generated di- 
rectly from the calibrated data stream, using 30-minute 
intervals of data for which the interferometer is operat- 
ing in what is known as science-mode. The choice of 30 



minutes is a tradeoff between intrinsic sensitivity, which 
increases with SFT length, and robustness against fre- 
quency drift during the SFT interval due to the Earth's 
motion, source spin-down, and non-stationarity of the 
data [7]. The requirement that each SFT contain contigu- 
ous data at nominal sensitivity introduces duty factor loss 
from edge effects, especially for the Livingston interfer- 
ometer (~20%) which had typically shorter contiguous- 
data stretches. In the end, the StackSlide and Hough 
searches used 1004 SFTs from HI and 899 from LI, 
the two interferometers with the best broadband sen- 
sitivty. For PowerFlux, the corresponding numbers of 
overlapped SFTs were 1925 and 1628. The Hough search 
also used 1063 H2 SFTs. In each case, modest require- 
ments were placed on data quality to avoid short periods 
with known electronic saturations, unmonitored calibra- 
tion strengths, and the periods immediately preceding 
loss of optical cavity resonance. 



B. Definitions And Notation 

Let N be the number of SFTs, T^„h the time-baseline 
of each SFT, and M the number of uniformly spaced 
data points in the time domain from which the SFT is 
constructed. If the time series is denoted by xj {j — 
0,1,2...M — 1), then our convention for the discrete 
Fourier transform is 



Xk — At Xje 

3=0 



-2mjk/M 



(12) 



where k = 0, 1, 2 . . . (M - 1), and At = T^^t^/M. For 
< k < M/2, the frequency index k corresponds to a 
physical frequency of fk = k/T^^^-,. 

In each method, the "power" (in units of spectral den- 
sity) associated with frequency bin k and SFT i is taken 
to be 



^ coh 



(13) 



It proves convenient to define a normalized power by 



Pk 



Pk 
SI 



(14) 



The quantity S], is the single-sided power spectral density 
of the detector noise at frequency fk, the estimation of 
which is described below. Furthermore, a threshold, Pthi 
can be used to define a binary count by jlOj : 



_ / 1 if Pi- > A> 
if p\ < A, 



(15) 



When searching for a signal using template A the de- 
tector antenna pattern and frequency of the signal are 
found at the midpoint time of the data used to gener- 
ate each SFT. Frequency dependent quantities are then 
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Quantity 



Description 



Pi Power for SFT i & template A 

Pi Normalized power for SFT i & template A 

rii Binary count for SFT i & template A 

Si Power spect. noise density for SFT i & template A 
F\. F+ at midpoint of SFT i for template A 

Fx at midpoint of SFT i for template A 

TABLE I: Summary of notation used. 



evaluated at a frequency index k corresponding to the bin 
nearest this frequency. To simplify the equations in the 
rest of this paper we drop the frequency index k and use 
the notation given in Table |T] to define various quantities 
for SFT i and template A. 



Basic StackSlide, Hough, and PowerFlux 
Formalism 



We call the detection statistics used in this search the 
"StackSlide Power" , P, the "Hough Number Count" , n, 
and the "PowerFlux Signal Estimator" , R. The basic 
definitions of these quantities are given below. 

Here the simple StackSlide method described in [15 
is used; the "StackSlide Power" for a given template is 
defined as 



1 



N-l 



i=0 



(16) 



This normalization results in values of P with a mean 
value of unity and, for Gaussian noise, a standard devi- 
ation of 1/^/N. Details about the value and statistics of 
P in the presence and absence of a signal are given in 
Appendix [b] and [TS]. 

In the Hough search, instead of summing the normal- 
ized power, the final statistic used in this paper is a 
weighted sum of the binary counts, giving the "Hough 
Number Count": 



N-l 
i=Q 

where the Hough weights are defined as 



(17) 



j{{nT + {K)'}, (18) 



binary count to have greater weight if the SFT i has a 
lower noise floor and if, in the time-interval correspond- 
ing to this SFT, the beam pattern functions are larger 
for a particular point in the sky. Note that the sensitiv- 
ity of the search is governed by the ratios of the different 
weights, not by the choice of overall scale. In the next sec- 
tion we show that these weights maximize the sensitivity, 
averaged over the orientation of the source. This choice 
of was originally derived in |16j using a different argu- 
ment and is similar to that used in the PowerFlux circular 
polarization projection described next. More about the 
Hough method is given in [7J [TD] . 

The PowerFlux method takes advantage of the fact 
that less weight should be given to times of greater noise 
variance or smaller detector antenna response to a signal. 
Noting that power estimated from the data divided by 
the antenna pattern increases the variance of the data at 
times of small detector response, the problem reduces to 
finding weights that minimize the variance, or in other 
words that maximize the signal-to-noise ratio. The re- 
sulting PowerFlux detection statistic is [17] , 



2 El-o'w,pj{F^y 



R = 



where the PowerFlux weights are defined as 



(20) 



(21) 



and where 



^pi\^2 _ ] (-P+)^ linear polarization ^^^^ 

^ y{F\Y + {FD'^ circular polarization 

As noted previously, the PowerFlux method searches us- 
ing four linear polarization projections and one circular 
polarization projection. For the linear polarization pro- 
jections, note that (F+)^ is evaluated at the angle -0, 
which is the same as {Fl^Y evaluated at the angle ip—Tr/A; 
for circular polarization, the value of (i^^)^-l-(F^)^ is in- 
dependent of 'if). Finally note that the factor of 2/T^„i, in 
Eq. ( 20 1 makes R dimensionless and is chosen to make it 



directly related to an estimate of the squared amplitude 
of the signal for the given polarization. Thus R is also 
called in this paper the "PowerFlux Signal Estimator". 
(See p/7j and Appendix [A] for further discussion.) 

We have shown in Eqs. (16l-(22l how to compute the 



detection statistic (or signal estimator) for a given tem- 
plate. The next section gives the details of the imple- 
mentation and pipelines used, where these quantities are 
calculated for a set of templates A and analyzed. 



and the weight normalization is chosen according to 



IMPLEMENTATIONS AND PIPELINES 



Af-l 
i=0 



(19) 



With this choice of normalization the Hough Number 
Count n lies within the range [0, A^]. Thus, we take a 



A. Running Median Noise Estimation 

The implementations of the StackSlide and Hough 
methods described below use a "running median" to esti- 
mate the mean power and, from this estimate, the power 
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spectral density of the noise, for every frequency bin of 
every SFT. PowerFlux uses a different noise decomposi- 
tion method described in its implementation section be- 
low. 

Note that for Gaussian noise, the single-sided power 
spectral density can be estimated using 



SI — 



(23) 



where the angle brackets represent an ensemble average. 
The estimation of must guard against any biases in- 
troduced by the presence of a possible signal and also 
against narrow spectral disturbances. For this reason the 
mean, is estimated via the median. We assume 

that the noise is stationary within a single SFT, but al- 
low for non-stationarities across different SFTs. In every 
SFT we calculate the "running median" of li^p for ev- 
ery 101 frequency bins centered on the k^^ bin, and then 
estimate [23l[24l|25] by dividing by the expected 

ratio of the median to the mean. 

Note, however, that in the StackSlide search, after 
the estimated mean power is used to compute S"^ in 



the denominator of Eq. ( 14 1 these terms are summed in 
Eq. 



( |16[ ) , while the Hough search applies a cutoff to ob- 
tain binary counts in Eq. ( 15 1 before summing. This re- 



sults in the use of a different correction to get the mean 
in the StackSlide search from that used in the Hough 
search. For a running median using 101 frequency bins, 
the effective ratio of the median to mean used in the 
StackSlide search was 0.691162 (which was chosen to nor- 
malize the data so that the mean value of the StackSlide 
Power equals one) compared with the expected ratio for 
an exponential distribution of 0.698073 used in the Hough 
search (which is explained in Appendix A of ^ ) . It is im- 
portant to realize that the results reported here are valid 
independent of the factor used, since any overall constant 
scaling of the data does not affect the selection of outliers 
or the reported upper limits, which are based on Monte 
Carlo injections subjected to the same normalization. 



B. The StackSlide Implementation 

1. Algorithm and parameter space 

The StackSlide method uses power averaging to gain 
sensitivity by decreasing the variance of the noise [121 HSl 
UM [T5] . Brady and Creighton [12 first described this ap- 
proach in the context of gravitational-wave detection as 
a part of a hierarchical search for periodic sources. Their 
method consists of averaging the power from a demodu- 
lated time series, but as an approximation did not include 
the beam pattern response of the detector. In Ref. [15] . 
a simple implementation is described that averages the 
normalized power given in Eq. ( 14 1 . Its extension to aver- 



aging the maximum likelihood statistic (known as the J-- 
statistic) which does include the beam pattern response 
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FIG. 3: Flow chart for the pipeline used to find the upper 
limits presented in this paper using the StackSlide method. 



is mentioned in Ref. [TS] (see also [6l[T0j[T8]), and further 
extensions of the StackSlide method are given in |T3] . 

As noted above, the simple StackSlide method given 
in [15 is used here and the detection statistic, called the 
"StackSlide Power" , is defined by Eq. ( 16 1 . The normal- 



ization is chosen so that the mean value of P is equal 
to 1 and its standard deviation is 1/y/N for Gaussian 
noise alone. For simplicity, the StackSlide Power signal- 
to-noise ratio (in general the value of P minus its mean 
value and then divided by the standard deviation of P) 
will be defined in this paper as (P— 1)^/N, even for non- 
Gaussian noise. 

The StackSlide code, which implements the method de- 
scribed above, is part of the C-based LSC Algorithms Li- 
brary Applications (LALapps) stored in the Iscsoft CVS 
repository f20J . The code is run in a pipeline with options 
set to produce the results from a search and from Monte 
Carlo simulations. Parallel jobs are run on computer 
clusters within the LSC, in the Condor environment [29], 
and the final post processing steps are performed using 
Matlab [30] • The specific StackSlide pipeline used to find 
the upper limits presented in this paper is shown in Fig. [3] 
The first three boxes on the left side of the pipeline can 
also be used to output candidates for follow-up searches. 

A separate search was run for each successive 0.25 Hz 
band within 50—1000 Hz. The spacing in frequency used 



is given by Eq. (111. The spacing in / was chosen as that 



which changes the frequency by one SFT frequency bin 
during the observation time Tobs, i.e., so that /Tobs = <5/. 
For simplicity T„b. = 2.778 x 10^ seconds ~ 32.15 days 
was chosen, which is greater than or equal to Tots for each 
interferometer. Thus, the / part of the parameter space 
was over-covered by choosing 



\Sf\ 



1 



2 X 10"^° Hz s 



(24) 
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Values of / in the range [—1 x 10~® Hz s~^,0 Hz s~^] 
were searched. This range corresponds to a search over 
51 values of /, which is the same as PowerFlux used in 
its low-frequency search (discussed in Section. VDl. 

The sky grid used is similar to that used for the all- 
sky search in [B], but with a spacing between sky-grid 
points appropriate for the StackSlide search. This grid is 
isotropic on the celestial sphere, with an angular spacing 
between points chosen for the 50-225 Hz band, such that 
the maximum change in Doppler shift from one sky grid 
point to the next would shift the frequency by half a bin. 
This is given by 



70 = 



O.bcSf „ „ . , /300Hz\ 

. ' = 9.3 X 10-3 , (25) 



145 146 147 14a 149 150 151 152 

Frequency (Hz) 



145 146 147 14a 149 150 151 152 153 154 155 

Frequency (Hz) 



where v is the magnitude of the velocity v of the detector 
in the SSB frame, and 9 is the angle between v and the 
unit-vector n giving the sky-position of the source. Equa- 
tions (24) and (25) are the same as Eqs. (19) and (22) in 
[7], which represent conservative choices that over-cover 
the parameter space. Thus, the parameter space used 
here corresponds to that in Ref. |3, adjusted to the S4 
observation time, and with the exception that a stere- 
ographic projection of the sky is not used. Rather an 
isotropic sky grid is used like the one used in [6j. 

One difficulty is that the computational cost of the 
search increases quadratically with frequency, due to the 
increasing number of points on the sky grid. To reduce 
the computational time, the sky grid spacing given in 
Eq. (25) was increased by a factor of 5 above 225 Hz. 



This represents a savings of a factor of 25 in computa- 
tional cost. It was shown through a series of simulations, 
comparing the upper limits in various frequency bands 
with and without the factor of 5 increase in grid spacing, 
that this changes the upper limits on average by less than 
than 0.3%, with a standard deviation of 2%. Thus, this 
factor of 5 increase was used to allow the searches in the 
225 — 1000 Hz band to complete in a reasonable amount 
of time. 

It is not surprising that the sky grid spacing can be in- 
creased, for at least three reasons. First, the value for Sdo 



given in Eq. ( 25 1 applies to only a small annular region on 



the sky, and is smaller than the average change. Second, 
only the net change in Doppler shift during the observa- 
tion time is important, which is less than the maximum 
Doppler shift due to the Earth's orbital motion during 
a one month run. (If the Doppler shift were constant 
during the entire observation time, one would not need 
to search sky positions even if the Doppler shift varied 
across the sky. A source frequency would be shifted by 
a constant amount during the observation, and would be 
detected, albeit in a frequency bin different from that 
at the SSB.) Third, because of correlations on the sky, 
one can detect a signal with negligible loss of SNR much 
farther from its sky location than the spacing above sug- 
gests. 



FIG. 4: The StackSlide Power for the 145-155 Hz band with 
no sliding. Harmonics of 1 Hz instrumental lines are clearly 
seen in HI (top) and LI (bottom). These lines are removed 
from the data by the StackSlide and Hough searches using the 
method described in the text, while PowerFlux search tracks 
these lines and avoids them when setting upper limits. 
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FIG. 5: The LI amplitude spectral density in a narrow fre- 
quency band estimated from 10 SFTs before and after the line 
cleaning used by the StackSlide pipeline. In the band shown, 
the 150 Hz bin, and one bin either side of this bin have been 
replaced with estimates of the noise based on neighboring 
bins. 



2. Line cleaning 

Coherent instrumental lines exist in the data which 
can mimic a continuous gravitational-wave signal for pa- 
rameter space points that correspond to little Doppler 
modulation. Very narrow instrumental lines are removed 
("cleaned") from the data. In the StackSlide search, a 
line is considered "narrow" if its full width is less than 
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IFO 


,/start 


,/stcp 


Num. 




A/right 


Description 




Hz 


Hz 




Hz 


Hz 




HI 


46.7 




1 


0.0 


0.0 


Cal. Line 


HI 


393.1 





1 


0.0 


0.0 


Cal. Line 


HI 


973.3 





1 


0.0 


0.0 


Cal. Line 


HI 


1144.3 





1 


0.0 


0.0 


Cal. Line 


HI 


0.0 


1.0 


1500 


0.0006 


0.0006 


1 Hz Comb 


LI 


54.7 




1 


0.0 


0.0 


Cal. Line 


LI 


396.7 




1 


0.0 


0.0 


Cal. Line 


LI 


1151.5 




1 


0.0 


0.0 


Cal. Line 


LI 


0.0 


1.0 


1500 


0.0006 


0.0006 


1 Hz Comb 



TABLE H: Instrumental lines cleaned during the StackSlide 
search. The frequencies cleaned are found by starting with 
that given in the first column, and then taking steps in fre- 
quency given in the second column, repeating this the num- 
ber of times shown in the third column; the fourth and fifth 
columns show how many additional Hz are cleaned to the 
immediate left and right of each line. 



Excluded Bands 


Description 


Hz 




[57, 63) 


Power lines 


[n60 - 1, n60 -f 1) n = 2 to 16 


Power line harmonics 


[340, 350) 


Violin modes 


[685, 690) 


Violin mode harmonics 


[693, 696) 


Violin mode harmonics 



TABLE HL Frequency bands excluded from the StackSlide 
search. 
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FIG. 6: Measure confidence vs. hg for an example band (140— 
140.25 Hz in HI). A best-fit straight line is used to find the 
value of ho corresponding to 95% confidence and to estimate 
the uncertainties in the results (see text). 



refer to resonant excitations of the steel wires that sup- 
port the interferometer mirrors.) These are shown in 
Table [Till While these bands can be covered by adjusting 
the parameters used to find outliers and set upper limits, 
we will wait for future runs to do this. 



3. Upper limits method 



5% of the 0.25 Hz band, or less than 0.0125 Hz. The line 
must also have been identified a priori as a known instru- 
ment artifact. Known lines with less than this width were 
cleaned by replacing the contents of bins corresponding 
to lines with random values generated by using the run- 
ning median to find the mean power using 101 bins from 
either side of the lines. This method is also used to esti- 
mate the noise, as described in Section [V A| 

It was found when characterizing the data that a comb 
of narrow 1 Hz harmonics existed in the HI and LI data, 
as shown in Fig.|4] Table shows the lines cleaned dur- 
ing the StackSlide search. As the table shows, only this 
comb of narrow 1 Hz harmonics and injected lines used 
for calibration were removed. As an example of the clean- 
ing process. Fig. [5] shows the amplitude spectral density 
estimated from 10 SFTs before and after line cleaning, 
for the band with the 1 Hz line at 150 Hz. 

The cleaning of very narrow lines has a negligible effect 
on the efficiency to detect signals. Very broad lines, on 
the other hand, cannot be handled in this way. Bands 
with very broad lines were searched without any line 
cleaning. There were also a number of highly disturbed 
bands, dominated either by the harmonics of 60 Hz power 
lines or by the violin modes of the suspended optics, that 
were excluded from the StackSlide results. (Violin modes 



After the lines are cleaned, the powers in the SFTs 
are normalized and the parameter space searched, with 
each template producing a value of the StackSlide Power, 
defined in Eq. (16 1. For this paper, only the "loudest" 
StackSlide Power is kept, resulting in a value Pmax for 
each 0.25 Hz band, and these are used to set upper lim- 
its on the gravitational-wave amplitude, Hq. (The loud- 
est coincident outliers are also identified, but none sur- 
vive as candidates after follow-up studies described in 



Sec. VII A 1 ) The upper limits are found by a series of 



Monte Carlo simulations, in which signals are injected 
in software with a fixed value for /ig, but with otherwise 
randomly chosen parameters, and the parameter space 
points that surround the injection are searched. The 
number of times the loudest StackSlide Power found dur- 
ing the Monte Carlo simulations is greater than or equal 
to Pmax is recorded, and this is repeated for a series of 
ho values. The 95% confidence upper limit is defined to 
be the value of ho that results in a detected StackSlide 
Power greater than or equal to Pmax 95% of the time. 
As shown in Fig. [3] the line cleaning described above 
is done after each injection is added to the input data, 
which folds any loss of detection efficiency due to line 
cleaning into the upper limits self-consistently. 

Figure [6] shows the measured confidence versus ho for 
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an example frequency band. The upper limit finding pro- 
cess involves first making an initial guess of its value, then 
refining this guess using a single set of injections to find 
an estimate of the upper limit, and finally using this es- 
timate to run several sets of injections to find the final 
value of the upper limit. These steps are now described 
in detail. 

To start the upper limit finding process, first an initial 
guess, hf^°^^, is used as the gravitational-wave amplitude. 
The initial guess need not be near the sought-after upper 
limit, just sufficiently large, as explained below. A single 
set of n injections is done (specifically n = 3000 was used) 
with random sky positions and isotropically distributed 
spin axes, but all with amplitude /iq"°''^. The output 
list of StackSlide Powers from this set of injections is 
sorted in ascending order and the O.OSn'th (specifically 
for n = 3000 the 150th) smallest value of the StackSlide 
Power is found, which we call Pq.os, Note that the goal is 
to find the value of Hq that makes Pq.os = ^maxi so that 
95% of the output powers are greater than the maximum 
power found during the search. This is what we call the 
95% confidence upper limit. Of course, in general Po.os 
will not equal Pmax unless our first guess was very lucky. 
However, as per the discussion concerning Eq. (B5), P—1 



is proportional to Hq (i.e, removing the mean value due 
to noise leaves on average the power due to the presence 
of a signal). Thus, an estimate of the 95% /ig confidence 
upper limits is given by the following rescaling of hf^'^''^, 



7 GSt 



(26) 



05 



Hq^^ , is found from a sin- 



Thus an estimated upper limit, 
gle set of injections with amplitude Hq^"^^; the only re- 
quirement is that Hq^^^^ is chosen loud enough to make 
-Po.05 > 1- 



It is found that using Eq. ( 26 1 results in a estimate of 



the upper limit that is typically within 10% of the final 
value. For example, the estimated upper limit found in 
this way is indicated by the circled point in Fig. |6] The 
value of hg^^ then becomes the first value for ho in a series 
of Monte Carlo simulations, each with 3000 injections, 
which use this value and 8 neighboring values, measur- 
ing the confidence each time. The Matlab ^30j polyfit 
and polyval functions are then used to find the best-fit 
straight line to determine the value of ho corresponding 
to 95% confidence and to estimate the uncertainties in 
the results. This is the final step of the pipeline shown 
in Fig. [3] 



C. The Hough Transform Implementation 

1. Description of Algorithm 

The Hough transform is a general method for pattern 
recognition, invented originally to analyze bubble cham- 
ber pictures from CERN [3TJ it has found many 



applications in the analysis of digital images ^33i . This 
method has already been used to analyze data from the 
second science run (S2) of the LIGO detectors [7] and a 
detailed description can be found in lOJ. Here we present 
only a brief description, emphasizing the differences be- 
tween the previous S2 search and the S4 search described 
here. 

The Hough search uses a weighted sum of the binary 
counts as its final statistic, as given by Eqs. ( 15 1 and ( 19 1 



In the standard Hough search as presented in |7[IT0j. the 
weights are all set to unity. The weighted Hough trans- 
form was originally discussed in |16j . The software for 
performing the Hough transform has been adapted to use 
arbitrary weights without any significant loss in compu- 
tational efficiency. Furthermore, the robustness of the 
Hough transform method in the presence of strong tran- 
sient disturbances is not compromised by using weights 
because each SFT contributes at most Wi (which is of 
order unity) to the final number count. 

The following statements can be proven using the 
methods of |10] . The mean number count in the ab- 
sence of a signal is n = Np, where N is the number 
of SFTs and p is the probability that the normalized 
power, of a given frequency bin and SFT defined by 
Eq. (14 1, exceeds a threshold pth, i-e., p is the proba- 



bility that a frequency bin is selected in the absence of 
a signal. For unity weighting, the standard deviation is 
simply cr = ^y Np{l — p). However, with more general 
weighting, it can be shown that a is given by 



a- v/||w||Ml-p), 



(27) 



where | |w| p = Si^o^ '^1- ^ threshold riji, on the number 
count corresponding to a false alarm rate is given by 



1,, ^Np+ V2||w||2p(l-P) erfc-i(2aH) 



(28) 



Therefore depends on the weights of the correspond- 
ing template A. In this case, the natural detection statis- 
tic is not the "Hough Number Count" n, but the signifi- 
cance of a number count, defined by 



(29) 



where n and a are the expected mean and standard de- 
viation for pure noise. Values of s can be compared di- 
rectly across different templates characterized by differ- 
ing weight distributions. 



The threshold p^^ (c.f. Eq. 15 1 is selected to give the 
minimum false dismissal probability (3^ for a given false 
alarm rate. In [7] it was shown that the optimal choice for 
Pth is 1.6 which correspond to a peak selection probability 
0.2. It can be shown that the optimal choice 



P - 



-pth 



is unchanged by the weights and hence p^-, — 1.6 is used 
once more [Ml. 

Consider a population of sources located at a given 
point in the sky, but having uniformly distributed spin 
axis directions. For a template that is perfectly matched 
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FIG. 7: The improvement in the significance as a function of 
the mismatch in the sky-position. A signal is injected in fake 
noise at a = 5 = and the weights are calculated at a = 5 = 
S9. The curve is the observed significance as a function of 59 
while the horizontal line is the observed significance when no 
weights are used. See main text for more details. 



in frequency, spin-down, and sky-position, and given the 
optimal peak selection threshold, it can be shown [M] 
that the weakest signal that can cross the threshold nth 
with a false dismissal probability /3h has an amplitude 



where 



ho = 3.38 5^/2 



S = erfc"^(2aH) 



1/2 



1 



X,, = 



herfc^^(2/3H). 

inf} ■ 



(30) 

(31) 
(32) 



As before, F'^ and i^* are the values of the beam pat- 
tern functions at the mid-point of the i*'' SFT. To derive 
( 30 1 wc have assumed that the number of SFTs N is 



sufficiently large and that the signal is weak [TD] . 

From (30) it is clear that the scaling of the weights 
does not matter; Wi — > kwi leaves ho unchanged for any 
constant k. More importantly, it is also clear that the 
sensitivity is best, i.e. hg is minimum, when w • X is 



maximum: 



cx X, . 



(33) 



This result is equivalent to Eq. (18 I. 

In addition to improving sensitivity in single- 
interferometer analysis, the weighted Hough method al- 
lows automatic optimal combination of Hough counts 
from multiple interferometers of differing senstivities. 

Ideally, to obtain the maximum increase in sensitiv- 
ity, we should calculate the weights for each sky-location 
separately. In practice, we break up the sky into smaller 
patches and calculate one weight for each sky-patch cen- 
ter. The gain from using the weights will be reduced if the 
sky patches are too large. From equation (32 1, it is clear 



only through the beam pattern functions. Therefore, the 
sky patch size is determined by the typical angular scale 
over which and i<"x vary; thus for a spherical detector 
using the beam pattern weights would not gain us any 
sensitivity. For the LIGO interferometers, we have in- 
vestigated this issue with Monte-Carlo simulations using 
random Gaussian noise. Signals are injected in this noise 
corresponding to the HI interferometer at a sky-location 
(q;o, (5o), while the weights are calculated at a mismatched 
sky-position {ao + S9, So + 66). The significance values are 
compared with the significance when no weights are used. 
An example of such a study is shown in Fig. [7] Here, we 
have injected a signal at a = (5 = 0, cosi — 0.5, zero 
spin-down, — ^ — 0, and a signal to noise ratio corre- 
sponding approximately to a 6-cr level without weights. 
The figure shows a gain of ~ 10% at 56' = 0, decreas- 
ing to zero at 69 w 0.3 rad. We get qualitatively similar 
results for other sky-locations, independent of frequency 
and other parameters. There is an additional gain due 
to the non-stationarity of the noise itself, which depends, 
however, on the quality of the data. In practice, we have 
chosen to break the sky up into 92 rectangular patches 
in which the average sky patch size is about 0.4 rad wide, 
corresponding to a maximum sky position mismatch of 
69 = 0.2 rad in Fig. [7] 



2. The Hough Pipeline 

The Hough analysis pipeline for the search and for 
setting upper limits follows roughly the same scheme as 
in [7 . In this section we present a short description of 
the pipeline, mostly emphasizing the differences from [7] 
and from the StackSlide and PowerFlux searches. As 
discussed in the previous subsection, the key differences 
from the S2 analysis pT are (i) using the beam-pattern 
and noise weights, and (ii) using SFTs from multiple in- 
terferometers. 

The total frequency range analyzed is 50-1000 Hz, with 
a resolution 6f — l/T^oh as in (111. The resolution in / is 
2.2 X 10~^° Hz s~^ given in (24 1, and the reference time 



for defining the spin-down is the start-time of the obser- 
vation. However, unlike StackSlide and PowerFluXj the 
Hough search is carried out over only 11 values of /, in- 
cluding zero, in the range [— 2.2x 10^^ Hz s^^, Hz s"-'^]. 
This choice is driven by the technical design of the cur- 
rent implementation, which uses look-up-tables and par- 
tial Hough maps as in [7]. This implementation of the 
Hough algorithm is efficient when analyzing all resolv- 
able points in /, as given in (24), but this approach is 



incompatible with the larger / step sizes used in the other 
search methods, which permit those searches to search a 
larger / range for comparable computational cost. 

The sky resolution is similar to that used by the Stack- 
Slide method for / < 225 Hz as given by (l25|). At fre- 



that the dependence of the weights on the sky-position is 



quencies higher than this, the StackSlide sky-resolution 
is 5 times coarser, thus the Hough search is analyzing 
about 25 more templates at a given frequency and spin- 
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TABLE IV: Instrumental lines cleaned during the Hough 
search that were not listed in Table [H] (see text). (*) These 
lines were removed only in the multi-interferometer search. 
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mean n = 202.7 and the standard deviation is given by 
Eq. (27 1. The standard deviation is computed from the 
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FIG. 8: Two example histograms of the normalized Hough 
number count compared to a Gaussian distribution for the HI 
detector in the frequency band 150-151 Hz. The upper figure 
corresponds to a a patch located at the north pole for the 
case in which the weights are used. The number of templates 
analyzed in this IHz band is of 11 x 10^, the number of SFTs 
1004, the corresponding mean n = 202.7 and cr = 12.94 is 
obtained from the weights. The lower figure corresponds to 
a patch at the equator using the same data. In this case the 
number of templates analyzed in this IHz band is of 10.5 x 10^, 
and its corresponding a — 14.96. 



down value. In each of the 92 sky patches, by means of 
the stereographic projection, the sky patch is mapped to 
a two dimensional plane with a uniform grid of that res- 
olution 69q. Sky Patches slightly overlap to avoid gaps 
among them (see [7j for further details) . 

Figure |8] shows examples of histograms of the number 
counts in two particular sky patches for the HI detector 
in the 150-151 Hz band. In all the bands free of instru- 
mental disturbances, the Hough number count distribu- 
tions follows the expected theoretical distribution, which 
can be approximated by a Gaussian distribution. Since 
the number of SFTs for HI is 1004, the corresponding 



weights w and varies among different sky patches because 
of varying antenna pattern functions. 

The upper limits on ho are derived from the loud- 
est event, registered over the entire sky and spin-down 
range in each 0.25 Hz band, not from the highest number 
count. As for the StackSlide method, we use a frequentist 
method, where upper limits refer to a hypothetical pop- 
ulation of isolated spinning neutron stars which are uni- 
formly distributed in the sky and have a spin-down rate 
/ imiformly distributed in the range [—2.2 x 10~^ Hz s~^, 
Hz s~^]. We also assume uniform distributions for the 
parameters cosi G [—1,1], "0 G [0,27r], and $0 G [0,27r]. 
The strategy for calculating the 95% upper limits is 
roughly the same scheme as in [7], except for the treat- 
ment of narrow instrumental lines. 

Known spectral disturbances are removed from the 
SFTs in the same way as for the StackSlide search. The 
known spectral lines are, of course, also consistently re- 
moved after each signal injection when performing the 
Monte-Carlo simulations to obtain the upper limits. 

The narrow instrumental lines "cleaned" from the SFT 
data are the same ones cleaned during the StackSlide 
search shown in Table [H] together with ones listed in Ta- 
ble |lVl The additional lines listed in Table ITVl are cleaned 
to prevent large artifacts in one instrument from increas- 
ing the false alarm rate of the Hough multi-interferometer 
search. Note that the LI 36.8725 Hz comb was eliminated 
mid-way through the S4 run by replacing a synthesized 
radio frequency oscillator for phase modulation with a 
crystal oscillator, and these lines were not removed in 
the Hough LI single-interferometer analysis. 

No frequency bands have been excluded from the 
Hough search, although the upper limits reported on the 
bands shown in Table [Till that are dominated by 60 Hz 
power line harmonics or violin modes of the suspended 
optics, did not always give satisfactory convergence to an 
upper limit. In a few of these very noisy bands, upper 
limits were set by extrapolation, instead of interpolation, 
of the Monte-Carlo simulations. Therefore the results 



14 



Noise decomposition 




Subtract medians 

r frequency, * 
: cumulate in TMedians 



Subtract medians 
accumulate in FMedians 



Coarse-grid SFT 
cutoff 4 
determination 



Loop over coarse grid 




Weighted means 




i 




Feldman-Cousins limits 






and statistics 








Update summary limits: 






maximum over skyband, 






frequency 








Log per-skyband information, 
optional output of 
sltymaps and plots 



FIG. 9: Flow chart for the pipeline used to find the upper 
limits presented in this paper using the PowerFlux method. 



reported on those bands have larger error bars. No pa- 
rameter tuning was performed on these disturbed bands 
to improve the upper hmits. 



D. The Powerflux Implementation 

The PowerFlux method is a variant on the StackSlide 
method in which the contributions from each SFT are 
weighted by the inverse square of the average spectral 
power density in each band and weighted according to 
the antenna pattern sensitivity of the interferometer for 
each point searched on the sky. This weighting scheme 
has two advantages: 1) variance on the signal strength 
estimator is minimized, improving signal-to-noise ratio; 
and 2) the estimator is itself a direct measure of source 
strain power, allowing direct parameter estimation and 
dramatically reducing dependence on Monte Carlo sim- 
ulations. Details of software usage and algorithms can 
be found in a technical document [17J. Figure |9] shows a 
flow chart of the algorithm, discussed in detail below. 



nominal 0.25 Hz band as a product of a spectral variation 
and a time variation across the data run. Specifically, 
for each 0.25 Hz band, a matrix of logarithms of power 
measurements across the 0.56 mHz SFT bins and across 
the SFT's of the run is created. Two vectors, denoted 
TMedians and FMedians, are initially set to zero and 
then iteratively updated according to the following 
algorithm: 

1. For each SFT (row in matrix), the median value 
(logarithm of power) is computed and then added 
to the corresponding element of TMedians while 
subtracted from each matrix element in that row. 

2. For each frequency bin (column in matrix), the me- 
dian value is computed and then added to the corre- 
sponding FMedians element, while subtracted from 
each matrix element in that column. 

3. The procedure repeats from step 1 until all medians 
computed in steps 1 and 2 are zero (or negligible). 

The above algorithm typically converges quickly. The 
size of the frequency band treated increases with central 
frequency, as neighboring bins are included to allow for 
maximum and minimum Doppler shifts to be searched in 
the next step. 

For stationary, Caussian noise and for noise that fol- 
lows the above assumptions of underlying factorized fre- 
quency and time dependence, the expected distribution 
of residual matrix values can be found from simulation. 
Figure [To] shows a sample expected residual power distri- 
bution following noise decomposition for simulated sta- 
tionary, Caussian data, along with a sample residual 
power distribution from the S4 data (0.25-Hz band of HI 
near 575 Hz, in this case) following noise decomposition. 
The agreement in shape between these two distributions 
is very good and is typical of the S4 data, despite some- 
times large variations in the corresponding TMedians and 
FMedians vectors, and despite, in this case, the presence 
of a moderately strong simulated pulsar signal (Pulsar2 
in Table |V]). 

The residuals are examined for outliers. If the largest 
residual value is found to lie above a threshold of 1.5, 
that corresponding 0.25 Hz band is flagged as containing 
a "wandering line" because a strong but drifting instru- 
mental line can lead to such outliers. The value 1.5 is 
determined empirically from Caussian simulations. An 
extremely strong pulsar could also be flagged in this way, 
and indeed the strongest injected pulsars are labelled as 
wandering lines. Hence in the search, the wandering lines 
are followed up, but no upper limits are quoted here for 
the affected bands. 



1. Noise decomposition 

Noise estimation is carried out through a 
time/frequency noise decomposition procedure in 
which the dominant variations are factorized within each 



2. Line flagging 

Sharp instrumental lines can prevent accurate noise es- 
timation for pulsars that have detected frequencies in the 
same 0.56 mHz bin as the line. In addition, strong lines 
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FIG. 10: Typical residual logarithmic power following noise 
decomposition for a sample 0.25-Hz band of HI data (crosses) 
near 575 Hz in a band containing an injected pulsar. The 
residual is defined as the difference between a measured power 
for a given frequency bin in a given 30-minute period and the 
value predicted by the FMedians and TMedians vectors. The 
smooth curve is for a simulation in Gaussian noise. 



PowerFlux algorithm computes a weighted sum of the 
strain powers, where the weighting takes into account 
the underlying time and spectral variation contained in 
TMedians and FMedians and the antenna pattern sensi- 
tivity for an assumed sky location and incident wave po- 
larization. Specifically, for an assumed polarization angle 
■0 and sky location, the following quantity is defined for 
each bin k of each SFT i: 



(34) 



where Fl is the T/i-dependent antenna pattern for the sky 



location, defined in Eq. (22 1. (See also Appendix [A| ) 
As in Sec. IV B| to simplify the notation we define 



— PJ{F^) as the value of Qi for SFT i and a given 
template A. 

For each individual SFT bin power measurement Pi, 
one expects an underlying exponential distribution, with 
a standard deviation equal to the mean, a statement that 
holds too for Qi. To minimize the variance of a signal 
estimator based on a sum of these powers, each contribu- 
tion is weighted by the inverse of the expected variance of 
the contribution. Specifically, we compute the following 
signal estimator: 



tend to degrade achievable sensitivity by adding excess 
apparent power in an affected search. In early LIGO sci- 
ence runs, including the S4 run, there have been sharp 
instrumental lines at multiples of 1 Hz or 0.25 Hz, arising 
from artifacts in the data acquisition electronics. 

To mitigate the most severe of these effects, the Pow- 
erFlux algorithm performs a simple line detection and 
flagging algorithm. For each 0.25 Hz band, the detected 
summed powers are ranked and an estimated Gaussian 
sigma computed from the difference in the 50% and 94% 
quantiles. Any bins with power greater than 5.0 a are 
marked for ignoring in subsequent processing. Specifi- 
cally, when carrying out a search for a pulsar of a nomi- 
nal true frequency, its contribution to the signal estima- 
tor is ignored when the detected frequency would lie in 
the same 0.56 mHz bin as a detected line. As discussed 
below, for certain frequencies, spin-downs and points in 
the sky, the fraction of time a putative pulsar has a de- 
tected frequency in a bin containing an instrumental line 
can be quite large, requiring care. The deliberate ignor- 
ing of contributing bins affected by sharp instrumental 
Unes does not lead to a bias in resulting limits, but it 
does degrade sensitivity, from loss of data. In any 0.25 
Hz band, no more than five bins may be flagged as lines. 
Any band with more than five line candidates is exam- 
ined manually. 



3. Signal estimator 

Once the noise decomposition is complete, with esti- 
mates of the spectral noise density for each SFT, the 



R 




(35) 
(36) 



where Pi and Qi are the expected uncorrected and 
antenna-corrected powers of SFT i averaged over fre- 
quency. Since the antenna factor is constant in this av- 
erage, — Pi/{F^)^- Furthermore, is a estimate of 
the power spectral den sity of the noise. The replacement 
Pi ^ S, gives Eq. 

Note that for an SFT i with low antenna pattern sen- 
sitivity \F^\, the signal estimator receives a small con- 
tribution. Similarly, SFT's i for which ambient noise 
is high receive small contributions. Because computa- 
tional time in the search grows linearly with the number 
of SFT's and because of large time variations in noise, 
it proves efficient to ignore SFT's with sky-dependent 
and polarization-dependent effective noise higher than a 
cutoff value. The cutoff procedure saves significant com- 
puting time, with negligible effect on search performance. 

Specifically, the cutoff is computed as follows. Let aj 
be the ordered estimated standard deviations in noise, 
taken to be the ordered means of Qi = jr^T,kQ%, where 
fcmax is the number of frequency bins use^ in the search 
template. Define jopt to be the index jmax for which 

1 /x;^»-cr2 is minimized. Only SFT's 



the quantity 

for which aj < "iaj^.^^^ are used for signal estimation. In 
words, Jopt defines the last SFT that improves rather 
than degrades signal estimator variance in an unweighted 
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FIG. 11: Sky map of run-summed PowerFlux weights for a 
0.25-Hz band near 575 Hz for one choice of hnear polarization 
in the S4 HI data. The normalization corresponds roughly to 
the effective number of median-noise SFT's contributing to 
the sum. 



mean. For the weighted mean used here, the effective 
noise contributions are allowed to be as high as twice the 
value found for jopt- The choice of 2(7^^^^ is determined 
empirically. 

The PowerFlux search sets strict, frequentist, all-sky 
95% confidence-level upper limits on the flux of gravita- 
tional radiation bathing the Earth. To be conservative in 
the strict limits, numerical corrections to the signal esti- 
mator are applied: 1) a factor of 1/ cos(7r/8) = 1.082 for 
maximum linear polarization mismatch, based on twice 
the maximum half-angle of mismatch (see Appendix |A]) 
and 2) a factor of 1.22 for bin-centered signal power loss 
due to Hann windowing (applied during SFT generation) ; 
and 3) a factor of 1.19 for drift of detected signal fre- 
quency across the width of the 0.56 mHz bins used in 
the SFT's. Note that the use of rectangular windowing 
would eliminate the need for correction 2) above, but 
would require a larger correction of 1.57 for 3) 

Antenna pattern and noise weighting in the PowerFlux 
method allows weaker sources to be detected in certain 
regions of the sky, where run-averaged antenna patterns 
discriminate in declination and diurnal noise variations 
discriminate in right ascension. Figure pT] illustrates the 
resulting variation in effective noise across the sky for a 
0.25-Hz HI band near 575 Hz for the circular polariza- 
tion projection. By separately examining SNR, one may 
hope to detect a signal in a sensitive region of the sky 
with a strain significantly lower than suggested by the 
strict worst-case all-sky frequentist limits presented here, 
as discussed below in section IVIDI Searches are carried 
out for four linear polarizations, ranging over polariza- 
tion angle from i/j = to tp = f"" in steps of tt/S and for 
(unique) circular polarization. 

A useful computational savings comes from defining 
two different sky resolutions. A "coarse" sky gridding is 
used for setting the cutoff value defined above, while fine 
grid points are used for both frequency and amplitude 
demodulation. A typical ratio of number of coarse grid 
points to number of fine grid points used for Doppler 
corrections is 25. 



4-. Sky banding 

Stationary and near-stationary instrumental spectral 
lines can be mistaken for a periodic source of gravita- 
tional radiation if the nominal source parameters are con- 
sistent with small variation in detected frequency during 
the time of observation. The variation in the frequency 
at the detector can be found by taking the time derivative 
of Eq. which gives, 

The detector's acceleration, a in this equation is domi- 
nated by the Earth's orbital acceleration aEarth, since the 
diurnal part of the detector's acceleration is small and 
approximately averages to zero during the observation. 
Thus, it should be emphasized that a single instrumental 
line can mimic sources with a range of slightly different 
frequencies and assumed different positions in the sky 
that lie in an annular band. For a source / assumed to 
be zero, the center of the band is defined by a circle 90 
degrees away from the direction of the average acceler- 
ation of the Earth during the run where aEarth • n = 0, 
i.e., toward the average direction of the Sun during the 
run. For source spin-downs different from zero, there can 
be a cancellation between assumed spin-down (or spinup) 
that is largely cancelled by the Earth's average accelera- 
tion, leading to a shift of the annular region of apparent 
Doppler stationarity toward (away from) the Sun. 

A figure of merit found to be useful for discriminating 
regions of "good" sky from "bad" sky (apparent detected 
frequency is highly stationary) is the "S* parameter" : 

S = /+[(J^ X VEarth/c) •n)]/o, (38) 

where is the Earth's angular velocity vector about the 
solar system barycenter. The term x VEarth is a mea- 
sure of the Earth's average acceleration during the run, 
where VEarth is taken to be the noise-weighted velocity 
of the HI detector during the run. Regions of sky with 
small |5| for a given / and / have stationary detected 
frequency. As discussed below in section |VID[ such re- 
gions are not only prone to high false-alarm rates, but 
the line fiagging procedure described in section |VD 2| 
leads to systematically underestimated signal strength 
and invalid upper limits. Hence limits are presented here 
for only sources with l^l greater than a threshold value 
denoted S'largo- The minimum acceptable value chosen 
for S large is fouud from software signal injections to be 
1.1 X 10~^ Hz s~^ for the 1-month S4 run and can be 
understood to be 

rr -^occupied bins /'Q0^ 

"Jlarge — Tf^ j [•J^) 

^ ohs ' ^ coh 

where iVoccupiod bins ~ 5 is the minimum total num- 
ber of 0.56 mHz detection bins occupied by the source 
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FIG. 12: S4 sky band regions (good - light gray, bad for LI 
- medium or dark gray, bad for HI & LI - dark gray) for a 
source frequency / = 100 Hz and three different assumed spin- 
down choices: a) zero; b) —3 x 10~^ Hz s~^; and c) — 1 x 10~* 
Hz s~^. The black circle indicates the average position of the 
Sun during the data run. 



Kiglirt Ascension 

FIG. 13: S4 sky band regions (good - light gray, bad for LI 
- medium or dark gray, bad for HI & LI - dark gray) for a 
source frequency / — 300 Hz and three different assumed spin- 
down choices: a) zero; b) —3 x 10~^ Hz s~^; and c) —1 x 10~* 
Hz s~^. The black circle indicates the average position of the 
Sun during the data run. 



during the data run for reliable detection. In prac- 
tice, we use still larger values for the HI interferome- 
ter {Siarge — 1-85 X 10^^ Hz s^^) and LI interferom- 
eter {Siarge = 3.08 x 10"^ Hz s^^) during the S4 run 
for the limits presented here because of a pervasive and 
strong comb of precise 1-Hz lines in both interferometers. 
These lines, caused by a GPS-second synchronized elec- 
tronic disturbance and worse in LI than in HI, lead to 
high false-alarm rates from that data for lower values of 
Siarge- For the frequency and spin-down ranges searched 
in this analysis, the average fractions of sky lost to the 
skyband veto are 15% for HI and 26% for LI. 

Figures [T2|T4| illustrate the variation in the fraction of 
sky marked as "bad" as assumed source frequency and 
spin-down are varied. Generally, at low frequencies, large 
sky regions are affected, but only for low spin-down mag- 
nitude, while at high frequencies, small sky regions are 
affected, but the effects are appreciable to larger spin- 
down magnitude. It should be noted that the annular 
regions of the sky affected depend upon the start time 
and duration of a data run. The longer the data run, the 
smaller is the region of sky for which Doppler stationar- 
ity is small. Future LIGO data runs of longer duration 
should have only small regions near the ecliptic poles for 
which stationary instrumental lines prove troublesome. 



5. Gnd-point upper limit determination 

An intermediate step in the PowerFlux analysis is the 
setting of upper limits on signal strength for each sky- 
point for each 0.56 mHz bin. The limits presented here 
for each interferometer are the highest of these interme- 
diate limits for each 0.25-Hz band over the entire "good" 
sky. The intermediate limits are set under the assump- 
tion of Gaussian residuals in noise. In brief, for each 
0.56 mHz bin and sky-point, a Feldman-Gousins |35j 
95% confidence-level is set for an assumed normal dis- 
tribution with a standard deviation determined robustly 
from quantiles of the entire 0.25 Hz band. The Feldman- 
Gousins approach provides the virtues of a well behaved 
upper limit even when background noise fluctuates well 
below its expectation value and of smooth transition be- 
tween I-sided and 2-sided limits, but in practice the high- 
est upper limit for any 0.25 Hz band is invariably the 
highest measured power plus 1.96 times the estimated 
standard deviation on the background power for that bin, 
corresponding to a conventional a priori I-sided 97.5% 
upper GL. A Kolmogorov-Smirnov (KS) statistic is com- 
puted to check the actual power against a Gaussian dis- 
tribution for each 0.25 Hz band. Those bands that fail the 
KS test value of 0.07 (> 5a deviation for the S4 data sam- 
ple) are flagged as "Non-Gaussian" , and no upper limits 
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FIG. 14: S4 sky band regions (good - light gray, bad for LI 

- medium or dark gray, bad for HI & LI - dark gray) for a 
source frequency / = 1000 Hz and three different assumed 
spin-down choices: a) zero; b) -3 X lO"'' Hz s-^ and c) 

— 1 X 10~* Hz s~^. The black circle indicates the average 
position of the Sun during the data run. 
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FIG. 15: Histogram of Feldman Cousins 95% confidence-level 
upper limits in a 0.25-Hz band near 149 Hz in S4 HI data. 
Each entry corresponds to the highest upper limit in the band 
for a single sky location. 



provided an end-to-end validation of the entire pipelines. 
The next subsections discuss the hardware injections, the 
validations of the three methods and their pipelines. The 
detection of the hardware injections also shows in dra- 
matic fashion that we can detect the extremely tiny sig- 
nals that the detectors were designed to find. 



on pulsars are quoted here for those bands, although a 
full search is carried out. Bands subject to violin modes 
and harmonics of the 60 Hz power mains tend to fail the 
KS test because of sharp spectral slope (and sometimes 
because non-stationarity of sharp features leads to poor 
noise factorization). 

Figure [15] provides an example of derived upper limits 
from one narrow band. The figure shows the distribu- 
tion of PowerFlux strain upper limits on linear polar- 
ization amplitude /iq™ for a sample 0.25 Hz band of S4 
HI data near 149 Hz. The highest upper limit found is 
3.35 X 10"^'' (corresponding to a worst-case pulsar up- 
per limit on /iq 

of 6.70 X 10-24). The bimodal distribu- 
tion arises from different regions of the sky with intrinsi- 
cally different antenna pattern sensitivities. The peak at 
2.8 X 10-24 corresponds to points near the celestial equa- 
tor where the run-averaged antenna pattern sensitivity is 
worst. 



VI. HARDWARE INJECTIONS AND 
VALIDATION 

All three methods discussed in this paper have under- 
gone extensive internal testing and review. Besides in- 
dividual unit tests of the software, hardware injections 



A. Hardware injections 

During a 15-day period in the S4 run, ten artificial iso- 
lated pulsar signals were injected into all three LIGO in- 
terferometers at a variety of frequencies and time deriva- 
tives of the frequency, sky locations, and strengths. Two 
additional artificial binary pulsar signals were injected 
for approximately one day. These hardware injections 
were implemented by modulating the interferometer mir- 
ror positions via signals sent to voice actuation coils sur- 
rounding magnets glued to the mirror edges. The in- 
jections provided an end-to-end validation of the search 
pipelines. Table [V] summarizes the nominal parameters 
used in the isolated-pulsar injections; the parameters are 
defined in section Hill 

Imperfect calibration knowledge at the time of these 
injections led to slightly different actual strain ampli- 
tude injections among the three LIGO interferometers. 
For the HI and LI comparisons between expected and 
detected signal strengths for these injections described 
in section |VIB| corrections must be applied for the dif- 
ferences from nominal amplitudes. The corrections are 
the ratios of the actuation function derived from final 
calibration to the actuation function assumed in the pre- 
liminary calibration used during the injections. For HI 
this ratio was independent of the injection frequency and 



19 



Name 


/o (Hz) 


df/dt (Hz s" 




a (radians) 


5 (radians) 


(radians) 


A+ 










PulsarO 


265.57693318 


-4.15 X 10" 


12 


1.248816734 


-0.981180225 


0.770087086 


4.0250 X 


10" 


25 


3.9212 X 10" 


25 


Pulsarl 


849.07086108 


-3.00 X 10" 


10 


0.652645832 


-0.514042406 


0.35603553 


2.5762 X 


10" 


24 


1.9667 X 10" 


24 


Pulsar2 


575.16356732 


-1.37 X 10" 


13 


3.75692884 


0.060108958 


-0.221788475 


7.4832 X 


10" 


24 


-7.4628 X 10 


-24 


PulsarS 


108.85715940 


-1.46 X 10" 


17 


3.113188712 


-0.583578803 


0.444280306 


1.6383 X 


10" 


23 


-2.6260 X 10 


-24 


Pulsar4 


1402.11049084 


-2.54 X 10" 


08 


4.886706854 


-0.217583646 


-0.647939117 


2.4564 X 


10" 


22 


1.2652 X 10" 


22 


PulsarS 


52.80832436 


-4.03 X 10" 


18 


5.281831296 


-1.463269033 


-0.363953188 


5.8898 X 


10" 


24 


4.4908 X 10" 


24 


Pulsar6 


148.44006451 


-6.73 X 10" 


09 


6.261385269 


-1.14184021 


0.470984879 


1.4172 X 


10" 


24 


-4.2565 X 10 


-25 


Pulsar? 


1220.93315655 


-1.12 X 10" 


09 


3.899512716 


-0.356930834 


0.512322887 


1.0372 X 


10" 


23 


9.9818 X 10" 


24 


PulsarS 


193.94977254 


-8.65 X 10" 


09 


6.132905166 


-0.583263151 


0.170470927 


1.5963 X 


10" 


23 


2.3466 X 10" 


24 


Pulsar9 


763.847316499 


-1.45 X 10" 


17 


3.471208243 


1.321032538 


-0.008560279 


5.6235 X 


10" 


24 


-5.0340 X 10 


-24 


PulsarlO 


501.23896714 


-7.03 X 10" 


16 


3.113188712 


-0.583578803 


0.444280306 


6.5532 X 10" 


23 


-1.0504 X 10 


-24 


Pulsarl 1 


376.070129771 


-4.2620 X 10 


-15 


6.132905166 


-0.583263151 


0.170470927 


2.6213 X 10" 


22 


-4.2016 X 10 


-23 



TABLE V: Nominal (intended) parameters for hardware injected signals, known as PulsarO to Pulsarll, for GPS reference time 
= 793130413 s (start of S4 run) at the SSB. These parameters are defined in section III As discussed in the text, imperfect 
calibration knowledge at the time of injections led to slightly different actual injected strain amplitudes among the three LIGO 
interferometers. The last two pulsars listed are binary system injections with additional orbital parameters not shown, which 
were injected during only the last day of the S4 run. 




equal to 1.12. For LI, this ratio varied slightly with fre- 
quency, with a ratio of 1.11 for all injected pulsars except 
Pulsarl (1.15) and Pulsar9 (1.18). 



B. StackSlide Validation 



Besides individual unit tests and review of each com- 
ponent of the StackSlide code, we have shown that sim- 
ulated signals are detected with the expected StackSlide 
Power, including the hardware injections listed in Ta- 
ble [V] Table VI shows the observed and injected SNR, 
and the square root of the observed and injected Stack- 
Slide Power, The percent difference of the latter is 
given, since this compares amplitudes, which are easier 
to compare with calibration errors. The observed val- 
ues were obtained by running the StackSlide code using 
a template that exactly matches the injection parame- 
ters, while the injected values were calculated using the 
parameters in Table [V] and the equations in Appendix 
[B] The SNR's of PulsarO, Pulsarl, PulsarS, and Pulsar6 
were too small to be detected, and Pulsar4 and Pulsar? 
were out of the frequency band of the all-sky search. Pul- 
sar2, PulsarS, and PulsarS were detected as outliers with 
SNR > ? (as discussed in Sec. 



RA (hours) 




RA (hours) 



FIG. 16: Detection of hardware injected Pulsar 2 by the 
StackSlide code in the HI (top) and LI (bottom) data. 



VII I 



while Pulsar9 was 

not loud enough to pass this requirement. In all cases 
the observed StackSlide Power agrees well with that pre- 
dicted, giving an end-to-end validation of the StackSlide 
code. 

As an example of an all-sky search for a band with an 
injection. Fig. [16] shows the detection of Pulsar2 for a 
search of the HI (top) and LI (bottom) data, and only 
during the times the hardware injections were running. 
Later, when the entire S4 data set was analyzed Pulsar2 
was still detected but with lower SNR, since this data in- 
cludes times when the hardware injections were absent. shown in Fig. 16 are 13.84 and 13.29. During the search 



Also note that, as explained in section |VD| because of 
strong correlations on the sky, a pulsar signal will be de- 
tected at many points that lie in an annular region in 
the sky that surrounds the point corresponding to the 
average orbital acceleration vector of the Earth, or its 
antipode. In fact, because of the large number of tem- 
plates searched, random noise usually causes the maxi- 
mum detected SNR to occur in a template other than 
the one which is closest to having the exact parameters 
of the signal. For example, for the exact template and 
times matching the Pulsar2 hardware injection, it was 
detected with SNR's of 8.92 and 8.20 in HI and LI, re- 
spectively, as given in Table |VI| while the largest SNR's 
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HI LI 
Observed Injected Observed Injected Percent Observed Injected Observed Injected Percent 



Pulsar 


SNR 


SNR 


VP 


VP 


Difference 


SNR 


SNR 


VP 


VP 


Difference 


PulsarO 


0.27 


0.23 


1.006 


1.005 


0.1% 


0.15 


0.13 


1.003 


1.003 


0.1% 


Pulsar 1 


1.62 


0.80 


1.035 


1.017 


1.7% 


0.27 


0.69 


1.006 


1.016 


-1.0% 


Pulsar2 


8.92 


8.67 


1.179 


1.175 


0.4% 


8.20 


9.34 


1.180 


1.203 


-1.9% 


PulsarS 


199.78 


174.72 


3.124 


2.943 


6.2% 


89.89 


104.76 


2.304 


2.454 


-6.1% 


Pulsar4 


2081.64 


1872.24 


9.607 


9.116 


5.4% 


1279.12 


1425.14 


7.895 


8.326 


-5.2% 


PulsarS 


0.05 


1.30 


1.001 


1.028 


-2.6% 


1.02 


0.44 


1.024 


1.010 


1.4% 


Pulsar6 


0.17 


2.94 


1.004 


1.063 


-5.5% 


2.90 


1.36 


1.067 


1.032 


3.4% 


Pulsar? 


6.25 


5.50 


1.129 


1.114 


1.3% 


6.07 


5.11 


1.136 


1.116 


1.8% 


PulsarS 


98.12 


96.21 


2.303 


2.285 


0.8% 


92.77 


103.45 


2.334 


2.441 


-4.4% 


Pulsar9 


6.68 


6.59 


1.137 


1.135 


0.2% 


2.61 


3.69 


1.061 


1.085 


-2.2% 



TABLE VI: Results of StackSlide analyses of the ten hardware injected continuous gravitational-wave signals from isolated 
neutron stars. 
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FIG. 17: Maximum significance as a function of frequency 
corresponding to the multi-interferometer search (using the 
data from the three detectors) and the HI and LI alone. 



of the full data set (including times when Pulsar2 was 
off) it was detected with SNR 11.09 and 10.71 in HI and 
LI, respectively. 
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FIG. 18: Maps of the Hough significance corresponding to 
the multi-interferometer case for Pulsar2, Pulsar3, Pulsar8 
and Pulsar9. The location of the injected pulsars are the cen- 
ters of the maps. For Pulsar2, PulsarS and Pulsar9, the maps 
correspond to the frequency and spin-down values closest to 
the real injected ones. For Pulsar8, we show the map con- 
taining the maximum significance value. The discrepancy in 
sky location is due to the mismatch in frequency and spin- 
down values between those used in the injections and those 
corresponding to the Hough map. 



C. Hough Validation 

Using the Hough search code, four hardware-injected 
signals have been clearly detected by analyzing the data 
from the interval when the injections took place. These 
correspond to Pulsar2, PulsarS, PulsarS and Pulsar9. For 
each of these injected signals, a small-area search (0.4 
radxO.4 rad) was performed, using a step size on the 
spin-down parameter of —4.2 x 10~^° Hz s~^. Given the 



large spin-down value of PulsarS (— S.65 x 10^^ Hz s^^), 
we have used 23 values of the spin-down spanning the 
range [—9.24 x 10~^ Hz s~-^, Hz s~-^] to search for this 
pulsar. Because of its large amplitude, PulsarS can be 
detected even with a large mismatch in the spin-down 
value, although at the cost of lower SNR. 

Figure [it] shows the significance maximized over differ- 
ent sky locations and spin-down values for the different 
frequencies. These four hardware injected pulsars have 
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Pulsar 


Detector 


fo range 


fo{max) 


Significance 






(Hz) 


(Hz) 




Pulsar2 


Multi-IFO 


575.15-575.18 


575.1689 


15.1195 




HI 


575.15-575.18 


575.1667 


11.1730 




LI 


575.15-575.18 


575.1650 


9.7635 


PulsarS 


Multi-IFO 


108.855-108.86 


108.8572 


39.1000 




HI 


108.855-108.86 


108.8572 


32.2274 




LI 


108.855-108.86 


108.8589 


19.2267 


PulsarS 


Multi-IFO 


193.932-193.945 


193.9411 


39.2865 




HI 


193.932-193.945 


193.9394 


27.9008 




LI 


193.932-193.945 


193.9400 


23.8270 


Pulsar9 


Multi-IFO 


763.83-763.87 


763.8511 


8.3159 




HI 


763.83-763.87 


763.8556 


6.1268 




LI 






5.4559 



TABLE VII: Results of the Hough search for the hardware 
injected signals for the multi-interferometer, HI and LI data. 
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Injected strain 



been clearly detected, with the exception of PulsarQ in 
the LI data. Pulsar9 is marginally visible using the HI 
data alone, with a maximum significance of 6.13, but 
when we combine the data from the three interferome- 
ters, the significance increases up to 8.32. Details are 
given in Table |VII[ including the frequency range of the 
detected signal, the frequency at which the maximum 
significance is obtained and its significance value. 



Figure 18 shows the Hough significance maps for the 
multi-interferometer case. The maps displayed corre- 
spond either to the frequency and spin-down values 
nearest to the injected ones, or to those in which the 
maximum significance was observed. The location of 
the injected pulsars correspond to the center of each 
map. Note that the true spin-down value of Pulsar8, 
—8.65 X 10^^ Hz s~^, lies between the parameter values 
-8.82 X 10~9 Hz s-i and -8.40 x 10"^ Hz s^^ of the 
nearest templates used. 



D. PowerFlux validation 

Several cross checks have been performed to validate 
the PowerFlux search algorithm. These validations range 
from simple and rapid Fourier-domain "power injections" 
to more precise time domain software simulations, to 
hardware signal injections carried out during data tak- 
ing. 

Signal strain power injections have been carried out 
as part of PowerFlux algorithm development and for pa- 
rameter tuning. These software injections involve super- 
imposing calculated powers for assumed signals upon the 
LIGO power measurements and carrying out searches. 
For computational speed, when testing signal detection 
efficiency, only a small region of the sky around the 
known source direction is searched. A critical issue is 
whether the strict frequentist limits set by the algorithm 
are sufficiently conservative to avoid undercoverage of 



FIG. 19: "Excess" (upper limit minus injected) strain plotted 
vs injected signal strain for sample PowerFlux HI elliptic- 
polarization near 140 Hz injections. 



the intended frequentist confidence band. We present 
here a set of figures that confirm overcoverage applies. 
Figure 19 shows the difference ("excess") between the 



Feldman-Cousins 95% confidence-level upper limit (con- 
ventional 97.5% upper limit) on strain and the injected 
strain for a sample of elliptic-polarization time-domain 
injections in the HI interferometer for the 140.50-140.75 
Hz band. Injection amplitudes were distributed logarith- 
mically, while frequencies, spin-downs, sky locations, and 
orientations were distributed uniformly. One sees that 
there is indeed no undercoverage (every excess strain 
value is above zero) over the range of injection ampli- 
tudes. Figure [20] shows the same "excess" plotted vs 
the injected spin-down value, where the search assumes 
a spin-down value of zero, and where the sample includes 
injections with actual spin-down values more than a step 
size away from the the assumed value for the search tem- 
plate. As one can see, in this frequency range, a spin- 
down stepsize of 1.0 x 10"^ Hz s~^ is safe (true spin- 
down no more 5.0 x 10"^*^ Hz away from the as- 



sumed search value). Figure 21 shows the "excess" plot 



ted vs the S parameter that discriminates between sky 
regions of low and high Doppler stationarity. As shown, 
a value of S'largc ~ 1 x 10^^ Hz s~^ is safe for these in- 
jections. For this search we have chosen 51 spin-down 
steps of 2 X 10"^° Hz s'^ for 50-225 Hz and 11 steps of 
1 X 10-9 Hz s~i for 200-1000 Hz. 

More computationally intensive full time-domain sig- 
nal injections were also carried out and the results found 
to be consistent with those from power injections, within 
statistical errors. 

In addition, the PowerFlux method was validated with 
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Detected /o (Hz) ho upper limit Det. polarization Detected SNR 



Pulsar 


/o (Hz) 


HI 


T 1 TniP hr. 




HI 




T 1 




HI 


T 1 


HI 


T 1 


Pulsar2 


575.164 


575.161 


575.164 8.04 x 10" 


24 


3.18 X 10" 


-23 


2.16 X 10" 


23 


circular 


circular 


16.59 


15.33 


PulsarS 


108.857 


108.858 


108.858 3.26 x 10" 


23 


3.92 X 10" 


-23 


3.36 X 10" 


23 


circular 


linear 


328.59 


209.99 


Pulsar4 


1402.110 


1402.111 


1402.113 4.56 X 10" 


22 


6.50 X 10" 


-22 


5.32 X 10" 


22 


linear 


circular 


2765.71 


1651.82 


Pulsar? 1220.933 


1220.933 


- 1.32 X 10" 


23 


3.56 X 10" 


-23 


2.88 X 10" 


23 


circular 




8.89 




PulsarS 


193.950 


193.951 


193.948 3.18 x 10" 


23 


4.18 X 10" 


-23 


3.52 X 10" 


23 


linear 


circular 


289.11 


292.13 


Pulsar9 


763.847 


763.849 


- 8.13 X 10" 


24 


1.69 X 10" 


-23 


1.97 X 10" 


23 


circular 




8.18 





TABLE VIII: Results of PowerFlux analysis of the six S4 hardware pulsar injections for which there is detection (SNR>7). 
Shown are the true nominal pulsar frequency at the start of the run (SSB frame), the frequency in each interferometer for 
detected signals, the true ho value of the injection, the worst-case upper limit from each interferometer, the polarization state 
for which the SNR is maximum in each interferometer, and the SNR of detected candidates. 




Spindown 




FIG. 20: "Excess" (upper limit minus injected) strain plotted 
vs injected signal spin-down for sample PowerFlux HI elliptic- 
polarization near 140 Hz injections. 



FIG. 21: "Excess" (upper limit minus injected) strain plotted 
vs S parameter defined in text, where values greater than 
8 X 10"^* have been "capped" at that ceiling value. 



VII. RESULTS 



the hardware signal injections summarized in Table |V] 
The PowerFlux algorithm was run on all 10 isolated pul- 
sars, including two outside the 50-1000 Hz search region, 
and results found to agree well with expectation for the 
strengths of the signals and the noise levels in their bands. 
Table |VIII| shows the results of the analysis for the six 
pulsars for which a detection with SNR>7 is obtained 
by PowerFlux for one or both of the 4-km interferome- 
ters. Figure [22] shows a sky map of PowerFlux -0 = 
polarization SNR for the 0.25 Hz band containing pulsar 
2 (575.16 Hz). 



All three methods described in Sections II VI and IVl have 
been applied in an all-sky search over a frequency range 
50-1000 Hz. As described below, no evidence for a grav- 
itational wave signal is observed in any of the searches, 
and upper limits on sources are determined. For the 
StackSlide and Hough methods, 95% confidence-level fre- 
quentist upper limits are placed on putative rotating 
neutron stars, assuming a uniform- sky and isotropic- 
orientation parent sample. Depending on the source loca- 
tion and inclination, these limits may overcover or under- 
cover the true 95% confidence-level band. For the Pow- 
erFlux method, strict frequentist upper limits are placed 
on linearly and circularly polarized periodic gravitational 
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FIG. 22: Sample sky map of Feldman Cousins upper limits 
on circularly polarized strain for a 0.25-Hz band containing 
hardware- injected Pulsar 2 at 575.16 Hz. Only the data (half 
the run) during which the pulsar injection was enabled has 
been analyzed for this plot. The injected pulsar {hg = 8.0 x 
10~^*) stands out clearly above background. (Right ascension 
increases positively toward the left and declination toward the 
top of the sky map.) 
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FIG. 23: The loudest observed StackShde Power for HI (top) 
and LI (bottom). Frequency bands with the harmonics of 60 
Hz and the violin modes have been removed. 



wave sources, assuming worst-case sky location, avoiding 
undercoverage. The limits on linear polarization are also 
re-interpreted as limits on rotating neutron stars, assum- 
ing worst-case sky location and worst-case star inclina- 
tion. The following subsections describe these results in 
detail. 
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FIG. 24: The loudest observed StackSlide Power for HI (top) 
and LI (bottom) with a simple veto applied: only outliers 
in each 0.25 Hz band with SNR > 7 in both interferometers 
that have a fractional frequency difference < 2.2 x 10~* are 
kept. These are shown against the background results that 
have SNR < 7 in both interferometers. Frequency bands with 
the harmonics of the 60 Hz and the violin modes have also 
been removed. 





fni (Hz) 


hi (Hz) 


HI SNR LI SNR 


Comment 


1 


78.618889 


78.618889 


14.82 


13.58 


Inst. 


Lines 


2 


108.856111 


108.856111 


152.11 


69.79 


HW Inj 


Pulsar3 


3 


193.947778 


193.949444 


121.89 


125.75 


HW Inj 


Pulsar8 


4 


244.148889 


244.157778 


9.00 


22.89 


Inst. 


Lines 


5 


375.793889 


375.806667 


11.68 


27.09 


HW Inj. 


Pulsar 11 


6 


376.271111 


376.281667 


7.47 


9.46 


HW Inj. 


Pulsar 11 


7 


575.162778 


575.153333 


11.09 


10.71 


HW Inj 


Pulsar2 


8 


575.250000 


575.371667 


7.49 


7.51 


Inst. 


Lines 


9 


575.250000 


575.153333 


7.49 


10.71 


Inst. & 


Pulsar2 


10 


580.682778 


580.734444 


7.02 


7.19 


Inst. 


Lines 


11 


912.307778 


912.271111 


7.02 


7.37 


Inst. 


Lines 


12 


988.919444 


988.960556 


9.56 


9.75 


Inst. 


Lines 


13 


988.919444 


989.000000 


9.56 


8.12 


Inst. 


Lines 


14 


993.356111 


993.523333 


7.08 


7.12 


Inst. 


Lines 



TABLE IX: StackSlide outliers with SNR > 7 in both inter- 
ferometers, with fraction difference in frequency less than or 
equal to 2.2 x 10~^, and after removal of the bands with 60 Hz 
harmonics and the violin modes. 



A. StackSlide Results 

1. Loudest powers and coincidence outliers 

The StackSlide method was applied to the S4 HI and 
LI data set, as given in Sec. VB As described in that 
section, only the loudest StackSlide Power was returned 
from a search of the entire sky, the range of the fre- 
quency's time derivative, [—1 x 10~®,0] Hz s^^, and for 



each 0.25 Hz band within 50 — 1000 Hz. The results are 
shown in Fig. [23] 

Many of the StackSlide results have power greater 
than expected due to random chance alone (for Gaus- 
sian noise). To identify the most interesting subset of 
these cases, a simple coincidence test was applied: only 
results with an SNR greater than 7 in both HI and LI 
and with a fractional difference in frequency, measured 
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FIG. 25: The StackSlide Power vs. frequency for HI (top), 
LI (middle) and H2 (bottom) using the sky position and the 
/ value of the template that gives the outlier in HI, for outlier 
number 2 given in Table|lX] Comparing with Tables|V]and|Vl] 
this outlier is identified as due to hardware injection PulsarS. 



FIG. 26: The StackSlide Power vs. frequency for HI (top), 
LI (middle) and H2 (bottom) using the sky position and the 
/ value of the template that gives the outlier in HI for outlier 
number 1 given in Table [iXl 



in the SSB, less than or equal to 2.2 x 10~^ were iden- 
tified as outliers for further follow-up. The requirement 
on frequency agreement comes from the worst-case sce- 
nario where a signal is detected on opposite sides of the 
sky with opposite Doppler shifts of 1 + v/c and 1 — w/c, 
giving a maximum fraction difference in the detected fre- 
quency at the SSB of 2v/c < 2.2 x IQ-*. The results 
after applying this simple coincidence test are shown in 
Fig. [24] The outliers that passed the test are shown in 
Tablelxl 

Note that the coincidence test used on the StackSlide 
results is very conservative in that it only covers the 
worst-case frequency difference, and makes no require- 
ment on consistency in sky position or the frequency's 
time derivative. However it is meant to find only the 
most prominent outliers. Since an automated follow-up 
of possible candidates is not yet in place, the follow-up is 
carried out manually. This dictated using a large thresh- 
old on SNR. Also, since the false dismissal rate of the 
coincidence test used was not determined (though it is 
assumed to be essentially zero) it is not used in this pa- 
per when setting upper limits. Monte Carlo studies will 
be needed to find appropriate thresholds on SNR and the 
size of coincidence windows, so that proper false alarm 
and false dismissal rates can be determined; such studies 
will be carried out when analyzing future data sets. 

Three types of qualitative follow-up tests were per- 
formed on each of the outliers in Table |IX| First, using 
the sky position and the / value of the template that 
gives the outlier in HI, the StackSlide Power was found 
using the same values for these in LI and H2 for a fre- 
quency band around that of the outlier in HI. For a 
fixed sky position and /, a true gravitational- wave sig- 



nal should show up in all three detectors as a narrow 
line at nearly the same frequency (though with an SNR 
corresponding to half the length displacement in H2 com- 
pared with that in HI and LI). Second, the StackSlide 
Power was computed for the frequency bands containing 
the outliers, with sliding turned off. If an instrumental 
line is the underlying cause of the outlier, a stronger and 
narrower peak will tend to show up in this case. Third, 
the StackSlide Power was computed for each HI outlier 
template, using half (and some other fractions) of the 
data. This should reduce the SNR of a true signal by 
roughly the square root of the fractional reduction of the 
data, but identify transient signals, which would fail this 
test by showing up in certain stretches of the data with 
more SNR while dissappearing in other stretches. This 
would be true of the hardware injections which were not 
always on during the run, or temporary disturbances of 
the instrument which appear to look like signals only for 
limited periods of time. (The search described here was 
not designed to find truely transient gravitational- wave 
signals.) 

The follow-up tests on the outliers given in Table |IX| 
found that none is qualitatively consistent with a true 
gravitational- wave signal. The three loudest hardware 
injections of periodic gravitational waves from fake iso- 
lated sources were found (indicated as PulsarS, PulsarS, 
and Pulsar2), as well as interference from a fake source in 
a binary system (Pulsarll). All of the outliers due to the 
hardware injections show up in the HI template as rela- 
tively narrow lines in all three detectors, for example as 
shown in Fig. 25 These outliers, on the other hand, fail 



the third test when looking at times the hardware injec- 
tions were turned off. In particular, this test, along with 
the frequencies in Table [V] confirms the identification of 
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Detector 


Band (Hz) 


,95% 


HI 
LI 


139.50-139.75 
140.75-141.00 


4.39 X 10"^* 
5.36 X 10-2* 



TABLE X: Best StackSlide all-sky ho upper limits obtained 
on the strength of gravitational waves from isolated neutron 
stars. 



outliers 5 and 6 as due to Pulsar 11. The other hardware 
injections also are identified as such via their detected 
frequencies in Table |V] and SNRs in Table |VI[ In com- 
parison, none of the other outliers qualitatively passes 



the first test, for example as shown in Fig. 26 The sec- 
ond test was less conclusive, since some of the outliers 
lie at points on the sky that receive little Doppler mod- 
ulation, but based on the first test we conclude that the 
remaining outliers are only consistent with instrumental 
line artifacts. These results are summarized in column 
six of Table [TX] In future searches, tests of the type used 
here should be studied using Monte Carlo simulations, to 
make them more quantitative. 



2. StackSlide upper limits 

The StackSlide 95% confidence upper limits on ho are 
shown as crosses for HI (top) and LI (bottom) respec- 
tively in Fig. [27] while the solid curves in this figure 
show the corresponding characteristic amplitudes given 
by Eq. (Bill in Appendix Ib| The characteristic ampli- 



tudes were calculated using an estimate of the noise from 
a typical time during the run, but include bands with the 
power line and violin line harmonics which were excluded 
from the StackSlide search. The best upper limits over 
the entire search band are given in Table |X] The uncer- 
tainties in the upper limits and confidence due to the 
method used are less than or equal to 3% and 5.3% re- 
spectively; random and systematic errors from the cali- 
bration increase these uncertainties to about 10%. 



B. Hough results 

1. Number Counts 

For the S4 data set, there are a total of = 2966 
SFTs from the three interferometers, giving an expected 
average number count for pure noise of n = Np ~ 593. 
The standard deviation a now depends on the sky-patch 
according to (27). For reference, if we had chosen unit 



weights, the standard deviation assuming pure Gaussian 
noise would have been ~ 22 for the multi-interferometer 
search. To compare number counts directly across differ- 
ent sky-patches, we employ the significance s of a number 
count defined in Eq. ( 29 1 . 



live contributions to the total Hough number count, and 
whether any of the interferometers should be excluded 
from the search, or if all of them should be included. 
For this purpose, for the moment let us ignore the beam 
pattern functions and consider just the noise weighting: 
Wi oc l/S.j^. The relative contribution of a particular in- 
terferometer, say /, is given by the ratio 



ri 



E 



N 



/ = H1, LI, H2. 



(40) 



The numerator is a sum of the weights for the / in- 
terferometer while the denominator is the sum of all the 
weights. This figurc-of-merit incorporates both the noise 
level of data from an interferometer, and also its duty 
cycle as determined by the number of SFTs available for 
that interferometer. Figure [28] shows the relative contri- 
butions from HI, LI, and H2 for the duration of the S4 
run. From the plot, we see that HI clearly contributes the 
most. H2 contributes least at low frequencies while LI 
contributes least at higher frequencies. Hence all three 
LIGO interferometers are included in this search. For 
comparison purposes and for coincidence analysis, we 
have also analyzed the data from HI and LI separately. 

Figure [29] shows the result of the Hough search using 
data from all three LIGO interferometers, either com- 
bined in a multi-interferometer search, or just for HI 
and LI data. This figure shows the loudest significance 
in every 0.25 Hz band, maximized over all sky-positions, 
frequencies and spin-downs for the three searches. Line 
cleaning was used as described before. In the bands 
in which there are no spectral disturbances the signifi- 
cance distribution agrees very well with the theoretical 
expected distribution as was shown in Fig. Js] 



2. Study of coincidence outliers 

There are many outliers from the Hough search with 
significance values higher than expected for Gaussian 
noise, as shown in Fig.j29] Many of the large outliers cor- 
respond to well known instrumental artifacts described 
earlier, such as the power mains harmonics or the violin 
modes. 

Note the relation between significance and false alarm 
which can be derived from equations ( 28 1 and ( 29 1 for 
Gaussian noise: 



an = 0.5erfc(s/\/2 



(41) 



Since the three interferometers have different noise 
floors and duty factors, we would like to know their rela- 



To identify interesting candidates, we consider only those 
that have a significance greater than 7 in the multi- 
interferometer search (the most sensitive one). This is the 
same threshold considered by the StackSlide and Pow- 
erFlux searches. For the Hough search, this threshold 
corresponds to a false alarm rate of 1.3 x 10"^^. With 
this threshold, we would expect about 6 candidates in 
a 100 Hz band around 1 kHz for Gaussian noise, since 
the number of templates analyzed in a 1 Hz band around 
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FIG. 27: The solid curve shows the characteristic amphtude given by Eq. (Bll I and crosses show the measured upper hmits 
on ho for the StackSlide search of the HI (top) and LI (bottom) data. 



1 kHz is about n = 4.4 x 10"'^°. If we would like to set 
a different threshold in order to select, say one event in 
a 1 Hz band, then we should increase the false alarm to 
an = l/n = 2.2 x 10"". 

In order to exclude spurious events due to instrumen- 
tal noise in just one detector, we pass these candidates 
through a simple coincidence test in both the HI and the 
LI data. Since the single detector search is less sensi- 
tive than the multi-interferometer one, we consider events 
from HI and LI with a significance greater than 6.6, cor- 
responding to a false alarm rate of 2.0 x 10"^^. The 
numbers of templates analyzed using the HI or LI data 
are the same as for the multi-interferometer search. 

The coincidence test applied first in frequency is sim- 
ilar to the one described for the StackSlide search, us- 
ing a coincidence frequency window as broad as the size 
of the maximum Doppler shift expected at a given fre- 
quency. Of the initial 3800 0.25-Hz bands investigated, 
276 yielded outliers in the multi-interferometer search 
with a significance higher than 7. Requiring those bands 







Hough 


significance 






Band (Hz) 


Multi-IFO 


HI LI 


Comment 


1 


78.602-78.631 


12.466 


12.023 10.953 


Inst. Lines 


2 


108.850-108.875 


29.006 


23.528 16.090 


Inj. Pulsar3 


3 


130.402-130.407 


7.146 


6.637 6.989 


? 


4 


193.92-193.96 


27.911 


17.327 20.890 


Inj. Pulsar8 


5 


575.15-575.23 


13.584 


9.620 10.097 


Inj. Pulsar2 


6 


721.45-721.50 


8.560 


6.821 13.647 LI Inst. Lines 


7 


988.80-988.95 


7.873 


8.322 7.475 


Inst. Lines 



TABLE XI; Hough outliers that have survived the coincidence 
analysis in frequency, excluding those related to 60 Hz har- 
monics and the violin modes. 



(or neighboring bands) to have outliers in HI higher 
than 6.6, reduced by half the number of surviving bands. 
These remaining bands were studied in detail and, after 
eliminating power line harmonics and the violin modes, 
27 candidates remained. Applying again the same co- 
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FIG. 28: Relative contributions of the three interferometers 
in the Hough muhi-interferometer search. The noise weights 
are calculated in 1 Hz bands. 
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FIG. 29: The measured loudest significance in each 0.25 Hz 
from the Hough search of the multi- interferometer (top), HI 
(middle) and Ll(bottom) data. 



incidence test with the LI data, we are left with only 
7 coincidence outliers that are listed on Table IXII and 
displayed in Fig. 30 

Except for the third outlier, the coincidence can be 
attributed to instrumental lines in the detectors or to 
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FIG. 30: Hough significance of the outliers that have survived 
the coincidence analysis without considering the bands con- 
taminated with 60 Hz harmonics or the violin modes. Points 
are plotted only for multi-interferometer templates with sig- 
nificance greater than 7 and for single-interferometer tem- 
plates with significance greater than 6.6. 



Detector s fo (Hz) df/dt (Hz s"^) q (rad) 5 (rad) 



Multi-IFO 7.146 130.4028 
HI 6.622 130.4039 
HI 6.637 130.4050 
LI 6.989 130.4067 



-1.745 X 10"^ 0.8798 
-1.334 X 10"^ 
-1.334 X 10"^ 
-1.963 X 10"^ 



-1.2385 
2.1889 0.7797 
2.0556 0.6115 
1.1690 -1.0104 



TABLE XII: Parameters of the candidate events with a sig- 
nificance greater than 6.6 in the multi-interferometer, HI and 
LI data searches around the Hough outlier number 3. The pa- 
rameters correspond to the significance, frequency and spin- 
down for the reference time of the beginning of S4, and sky 
locations. 



the hardware pulsar injections. Table XII summarizes 
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Detector 


Band (Hz) 


.95% 




H1+H2+L1 


140.00-140.25 


4.28 X 


10" 


-24 


HI 


129.00-129.25 


5.02 X 


10" 


-24 


LI 


140.25-140.50 


5.89 X 


10" 


-24 



TABLE XIII: Best Hough all-sky upper limits obtained on the 
strength of gravitational waves from isolated neutron stars. 



the parameters of the third coincidence candidate in the 
130.40-130.41 Hz frequency band, including all the events 
that in any of the searches had a significance larger than 
6.6. As can be seen from the Table, the events from the 
different data sets correspond to widely separated sky 
locations. Hence no detections were made in the Hough 
search of the S4 data. 

In future searches we plan to use lower thresholds in the 
semi-coherent step in order to point to interesting areas 
in parameter space to be followed up, using a hierarchi- 
cal scheme with alternating coherent and semicoherent 
steps. In what follows we will concentrate on setting up- 
per limits on the amplitude ho in each of the 0.25 Hz 
bands. 



3. Upper limits 

As in the previous S2 Hough search [7], we set a pop- 
ulation based frequentist upper limit using Monte Carlo 
signal software injections. We draw attention to two im- 
portant differences from that analysis: 

• In [7] , known spectral disturbances were handled by 
simply avoiding all the frequency bins which could 
have been affected by Doppler broadening. Thus, 
the loudest event was obtained by excluding such 
frequency bins, and the subsequent Monte Carlo 
simulations also did not perform any signal injec- 
tions in these bins. Here we follow the same ap- 
proach as used in the StackSlide search; we use the 
spectral line removal procedure described in sec- 
tion |VB 1[ For consistency, the same line removal 
procedure is followed in the Monte Carlo simulation 
after every software injection. 

• Recall that the calculation of the weights depends 
on the sky-patch, and the search has been carried 
out by breaking up the sky in 92 patches. Thus, 
for every randomly injected signal, we calculate the 
weights corresponding to the center of the corre- 
sponding sky patch. The analysis of [7] did not use 
any weights and this extra step was not required. 

The 95% confidence all-sky upper limit results on Hq from 
the Hough search for the multi-interferometer, HI and 

These upper limits have 



LI data are shown in Fig. 31 



been obtained by means of Monte-Carlo injections in each 
0.25 Hz band in the same way as described in |7j. The 
best upper limit over the entire search band corresponds 



to 4.28 X 10 '^^ for the multi-interferometer case in the 
140.00 — 140.25 Hz band. The results are summarized in 
Table IXlIIl 

Let us now understand some features of the upper-limit 
results. First, it turns out that it is possible to accurately 
estimate the upper limits without extensive Monte Carlo 
simulations. From (30 1, and setting Wi oc Xi, we expect 



that the upper limits are: 



,95% 



IXI 



1/2 



(42) 



Recall that Xi contains contributions both from the sky- 
location-dependent antenna pattern functions and from 
the sky-location-independent noise floor estimates. How- 
ever, since we are setting upper limits for a population 
uniformly distributed in the sky, we might expect that 
the 5. are more important for estimating the value of 



hf^°. From Eq. (1321 and averaging over the sky we get 



|X||a 



N-1 

E 

j=0 



(43) 



and thus, up to a constant factor C, the estimated upper 
limits are given by 



,95% 



c 



1 



1/4 



s 



(44) 



The value of S is calculated from Eq. (31 1 using the false 
alarm an corresponding to the significance of the ob- 
served loudest event in a particular frequency band. The 
value of the false dismissal rate /3h corresponds to the 
desired confidence level of the upper limit (in this case 
95%). To show that such a fit is viable. Fig. 32 plots the 



value of the constant C appearing in the above equation 
for every 0.25 Hz frequency band, using the measured up- 
per limits. It turns out that C — 9.2 ± 0.5. The exact 
value of C depends on the interferometer and the search 
performed, but it is still found to lie within this range. 
This scale factor C = 9.2 ± 0.5 is about two times worse 
than we would expect if we were performing a targeted 
(multi-interferometer with weights) search with no mis- 
match. This factor of two is also in very good agreement 
with what was reported in the S2 search [7]. 

The utility of this fit is that having determined the 
value of C in a small frequency range, it can be extrap- 
olated to cover the full bandwidth without performing 
any further Monte Carlo simulations. Figure [33] plots 
the ratio of the measured upper limits to the estimated 
values showing the accuracy of the fit. The scale factors 
C used are 9.2 for the multi-interferometer search, 9.7 for 
HI and 9.3 for LI. The scale factors have been obtained 
in all cases by comparing the measured upper limits by 
means of Monte Carlo injections to the quantity h^'^° /C 
as defined in Equation (44), using the full bandwidth of 
the search. These estimated upper limits have an error 
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FIG. 31: The 95% confidence all-sky upper limits on ho from the Hough search of the multi-interferometer (top), HI (middle) 
and LI (bottom) data. These upper limits have been obtained by means of Monte-Carlo injections in each 0.25 Hz band. 



smaller than 5% for bands free of large instrumental dis- 
turbances. 

We conclude this section by quantifying the improve- 
ment in sensitivity caused by using the weights. Fig- 
ure 34 shows the comparison between the weighted and 



un-weighted results in the 800-900 Hz frequency range. 
The average improvement is ~9% in this band. It is 
easy to see that the improvement as compared to the un- 
weighted Hough search will be larger if the variation of 
iSj and the beam pattern functions is large across the 
SFTs. Since the variation in S.^ is larger in a multi- 
interferometer search, we expect this improvement to be 
much more significant in a multi-interferometer search. 
For the case of analyzing data from a single interferome- 
ter, for example HI, the improvement in the upper limits 
due to the weights turns out to be only ^ 6%. Also, the 
improvement can be increased by choosing smaller sky- 
patches so that the weight calculation is more optimal. 
In particular, if there would not be any sky mismatch in 
computing the weights, only due to the amplitude mod- 



ulation, i.e., in the presence of Gaussian and stationary 
noise, we would expect an average increase of sensitiv- 
ity of ~10%, and it could be up to ~12% for optimally 
oriented pulsars. These results have been verified exper- 
imentally by means of a set of Monte-Carlo tests [33] • 



C. PowerFlux results 

1. Single-interferometer results 

The PowerFlux method has been applied to the S4 
data sample in the range 50-1000 Hz. Five polarization 
projections are sampled for each grid point: four linear 
polarizations with ip = 0, 7r/8, 7r/4, 3 7r/8; and circu- 
lar polarization. For each sky grid point in the "good 
sky" defined above and each of the 501 frequency bins 
(there is slight overlap of 0.25 Hz bands), the Feldman- 
Cousins [3S] 95% CL upper limit is computed, as de- 



scribed in section VD5 for each polarization projec- 
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FIG. 32: Ratio of the upper limits measured by means 
of Monte-Carlo injections in the multi-interferometer H ough 
search to the quantity /tp^^'/C as defined in Equation (44 1. 



The value of S in equation ( 44 1 is computed using the false 
alarm an corresponding to the observed loudest event, in a 
given frequency band, and for a false dismissal rate /3h ~ 0.05, 
in correspondence to the desired confidence level of the upper 
limit. The comparison is performed in each 0.25 Hz band. 
Analysis of the full bandwidth, and also in different 100 Hz 
bands, yield a scale factor C to be 9.2 ± 0.5. 



FIG. 34: Comparison of the upper limits obtained using 500 
Monte Carlo injections with and without weights in 0.5 Hz 
bands for the Hough multi-interferometer search. The use of 
the weights improves the upper limits by a ~9% factor. 
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FIG. 33: Ratio of the 95% confidence all-sky upper limits 
on ho obtained from the Hough search by means of Monte 
Carlo injections to those predicted by Eq. ( |44[ | of the multi- 
interferometer (top), HI (middle) and Ll(bottom) data. The 
comparison is performed in 0.25 Hz bands. The scale factors 
C used are 9.2 for the multi-interferometer search, 9.7 for HI 
and 9.3 for LI. 



tion. Worst-case upper limits on linear polarization for 
each grid point and frequency are taken to be the high- 
est linear-polarization-projection strain limit divided by 
cos(7r/8) to correct for worst-case polarization mismatch. 
The highest limit for all frequency bins in the 0.25 Hz 
band and over all sampled sky points is taken to be the 
broad-sky limit for that 0.25 Hz band. Figures 35]|36 
show the resulting broad-sky limits on linearly polarized 
periodic sources from HI and LI. Bands flagged as non- 
Gaussian (instrumental artifacts leading to failure of the 
KS test) or near 60-Hz harmonics are indicated by color. 
The derived upper limits for these bands are considered 
unreliable. Diamonds indicate bands for which wander- 
ing instrumental lines (or very strong injected signals) 
lead to degraded upper limits. An exceedingly strong 
pulsar can be identified as a wandering line, and several 
strong hardware-injected pulsars are marked in the fig- 
ures as such. 



These limits on linearly polarized radiation and the 
corresponding limits on circularly polarized radiation 
can be interpreted as worst-case and best-case limits 
on a triaxial-ellipsoid, non-precessing neutron star, re- 
spectively, as discussed in Appendix |A] Multiplying the 
linear-polarization limits by a factor of two leads to the 
worst-case HI limits on ho shown in Figs. |37|^38| The 
circular-polarization limits require no scale correction. 
Note that the StackSlide and Hough HI limits shown 
on the same figure apply to a uniform-sky, uniform- 
orientation population of pulsars. 
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FIG. 35: PowerFlux limits on linearly polarized CW radiation 
amplitude for the HI data from the S4 run. Bands flagged 
as non-Gaussian (instrumental artifacts) or near 60-Hz har- 
monics, and for which derived upper limits are unreliable, are 
indicated by color. Diamonds indicate bands for which wan- 
dering instrumental lines (or very strong injected signals) lead 
to degraded upper limits. 
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FIG. 36: PowerFlux limits on linearly polarized CW radiation 
amplitude for the LI data from the S4 run, with the same 
color coding as in the preceding figure. 



2. Coincidence followup of loud candidates 

All outliers (SNR>7, diamonds, and non-Gaussian 
bands) in the single-interferometer analysis are checked 
for coincidence between HI and LI. In this followup, 
agreement is required in frequency to within 10 mHz, in 
spin-down to within 1 x 10^^° Hz s^^, and in both right 
ascension and declination to within 0.5 radians. The 
only surviving candidates are associated with hardware- 
injected pulsars 2, 3, 4, and 8 (see Table VHI), 1-Hz 



harmonics, violin modes, and instrumental lines in both 
detectors near 78.6 Hz (also seen in the StackSlide and 
Hough searches). The source of these lines remains un- 
known, but followup consistency checks described in sec- 
tion 



VH A rule out an astrophysical explanation. 



From this coincidence analysis, we see no evidence of 
a strong pulsar signal in the S4 data. It should be noted, 
however, that the SNR threshold of 7 is relatively high. 
A lower threshold and a more refined algorithm for loca- 
tion and frequency coincidence is under development for 
future searches. 



VIII. COMPARISON OF THE THREE 
METHODS 



Figures |37] and |38] show superimposed the final upper 
limits on ho from the StackSlide, Hough, and PowerFlux 
methods when applied to the S4 single-interferometer HI 
and LI data, respectively. As one might have expected, 
we see that the StackSlide and Hough population-based 
limits lie between the best-case and worst-case ho strict 
limits from PowerFlux. As indicated in Figs. 37 38] the 
Hough search sensitivity improves with the summing of 
powers from two or more interferometers. 

To be more precise as to expectations, we have di- 
rectly compared detection efficiencies of the three meth- 
ods in frequency bands with different noise character- 
istics. As discussed above, we expect overall improved 
performance of Powerfiux with respect to StackSlide and 
Hough, except possibly for frequency bands marked by 
extreme non-Gaussianity or non-stationarity, where the 
Hough integer truncation of extreme power outliers can 
provide more robustness. We do not consider computa- 
tional efficiency, which could play an important role in 
deciding which algorithm to use in computationally lim- 
ited hierarchical searches. 



A comparison is shown in Figs. 39 and 40 among the ef- 
ficiencies of the three methods for two particular 0.25 Hz 
bands for HI: 140.5-140.75 Hz and 357-357.25 Hz. The 
horizontal axis in each case is the /iq of Monte Carlo soft- 
ware injections with random sky-locations, spin-downs 
and orientations. The noise in the two bands have qual- 
itatively different features. The 140.5-140.75 Hz band is 
a typical "clean" band with Gaussian noise and no ob- 
servable spectral features. 



39 shows 



As expected. Fig. 
that the efficiency for the PowerFlux method is higher 
than that for StackSlide, while that of StackSlide is bet- 
ter than that for Hough. In other bands, where there are 
stationary spectral disturbances, we find that PowerFlux 
remains the most efficient method. 

The noise in the band 357-357.25 Hz is non-Gaussian 
and displays a large transient spectral disturbance, in 
addition to stationary line noise at 357 Hz itself. The 
stationary 357 Hz line was removed during the Stack- 
Slide and Hough searches, avoided during the PowerFlux 
search, and handled self-consistently during Monte Carlo 
software injections. In this band, the Hough transform 



32 



1fr-22 - 



5e-23 - 



26-23 - 



1fr-23 - 



5fr-24 - 




□ PowerFlux worst-case (linear pol) 

□ Hough (H1 onl/) 

■ StackSlide 

■ Hough (multi-IFO) 

■ PowerFlux best-case [circ pol) 



16-24 - 



200 



400 



eoo 



800 



1000 



Frequency (Hz) 



FIG. 37: HI Upper limits (95% CL) on ho from the three methods. The StackShde and Hough limits are population-based, 
while those from PowerFlux are strict and apply, respectively, to the most favorable and least favorable pulsar inclinations. 
Also shown are the multi-interferometer limits from the Hough search. 



method proves to be robust against transient noise, and 
more sensitive than the StackShde or PowerFlux imple- 
mentations (see Fig. 39 1. In fact, no PowerFlux upper 
limit is quoted for this band because of the large non- 
Gaussianity detected during noise decomposition. Note 
that the SNR thresholds used for StacksHde, Hough and 
PowerFlux in Fig. |40]are set to 6.3, 5.2 and 30, respec- 
tively, to match their loudest events in this band of the 
data. 



IX. SUMMARY, ASTROPHYSICAL REACH, 
AND OUTLOOK 



In summary, we have set upper limits on the strength 
of continuous-wave gravitational radiation over a range 
in frequencies from 50 Hz to 1000 Hz, using three differ- 
ent semi-coherent methods for summing of strain power 
from the LIGO interferometers. Upper limits have been 
derived using both a population-based method applica- 
ble to the entire sky and a strict method applicable to 
regions of the sky for which received frequencies were not 



stationary during the S4 data run. 

The limits have been interpreted in terms of ampli- 
tudes ho for pulsars and in terms of linear and circu- 
lar polarization amplitudes, corresponding to least fa- 
vorable and most favorable pulsar inclinations, respec- 
tively. As a reminder, sets of known instrumental spec- 
tral lines have been cleaned from the data prior to setting 
the population-based StackShde and Hough upper limits 
( Tables [h] [ml and|lV|, while regions of the sky (defined 
by cutoff values on the S parameter (Equations |38] and 
|39[ ) have been excluded in the strict PowerFlux upper 
limits. The numerical values of the upper limits can be 
obtained separately |36j. 

We have reached an important milestone on the road 
to astrophysically interesting all-sky results: Our best 
upper limits on ho are comparable to the value of a few 
times lO^^'' at which one might optimistically expect to 
see the strongest signal from a previously unknown neu- 
tron star according to a generic argument originally made 
by Blandford (unpublished) and extended in our previ- 
ous search for such objects in S2 data !6^. The value from 
Blandford's argument does not depend on the distance to 
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FIG. 38: LI Upper limits (95% CL) on ho from the three methods. The StackSlide and Hough limits are population-based, 
while those from PowerFlux are strict and apply, respectively, to the most favorable and least favorable pulsar inclinations. 
Also shown are the multi-interferometer limits from the Hough search. 
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FIG. 39: Comparison of StackSlide, Hough, and PowerFlux 
efficiencies (SNR > 7) vs injected strain amplitude ho for the 
band 140.50-140.75 Hz for HI. From left to right, the curves 
correspond to PowerFlux, StackSlide, and Hough. This band 
is typical of those without large outliers. 



the star or its ellipticity, both of which are highly uncer- 
tain. 

We find the next milestone by considering the maxi- 
mum distance to which a signal could be detected and 
the ellipticity needed to generate a signal of the required 
strength at that distance. Both quantities are of inter- 
est since there are theoretical limits on the ellipticity, 
and both quantities are functions of the gravitatipnal- 
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wave frequency / and its derivative /. Figure 
contour plot of both quantities simultaneously, which we 
explain here in more detail. The Hough transform multi- 
interferometer upper limits on ho are used for illustration 
because they fall in the middle of the range of values for 
the different searches (see Fig. 37 1. The maximum dis- 
tance d{f, f) is obtained by equating the 95% confidence 
upper limits on ho for the multiple-interferometer plot 
in Fig. |3T]to the spin-down limit given in Eq. (0. This 
tacitly assumes that / is entirely due to emission of grav- 
itational radiation, which implies the ellipticity given in 
Eq. (H regardless of the data and the distance to the 
source. If we relaxed this assumption, knowing that neu- 
tron stars spin down due to electromagnetic wave emis- 
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FIG. 40: Detection efficiency curves for the frequency band 
357-357.25 Hz, for HI. This band has a transient spectral dis- 
turbance affecting some of the SFTs. The Hough transform 
method proves to be robust against such non-stationarities 
and is more sensitive than StackSlide or PowerFlux in this 
band. The SNR thresholds used to generate these curves 
were 6.3, 5.2, and 30, respectively, for the StackSlide, Hough, 
and PowerFlux methods, where the StackSlide and Power- 
Flux thresholds correspond to the loudest candidates in that 
band in the data. 



sion, relativistic particle winds, and other factors as well, 
the maximum distance and required ellipticity for a given 
/ and / would both be reduced. The degree of reduction 
would, however, be highly uncertain. 

We can use the combined contour plot in Fig. |4l] to 
answer questions about the astrophysical significance of 
our results. Here we ask and answer several salient ques- 
tions. First, what is the maximum range of the Hough 
transform search? The answer is obtained from looking 
at the top of Fig. [41] We could detect isolated pulsars 
to about 1 kpc, but only for a star radiating at a fre- 
quency near 100 Hz and then only if that star has an 
ellipticity somewhat more than 10"'', which is allowed 
only in the most extreme equations of state [371 ISHJ 155] . 
Second, what is the maximum range of detection for a 
normal neutron star? Normal neutron stars are expected 
to have e < 10~^ based on theoretical predictions [40] . 
By tracing the e = 10"^ contour, we find that the max- 
imum range is about 50 pc at the highest frequencies 
(1 kHz), falling with frequency to less than 2 pc below 
100 Hz. Third, what is the maximum range for a recy- 
cled millisecond pulsar? Based on the observed sample 
[22] , recycled pulsars usually have small |/| values, cor- 
responding to tsd usually less than 10"^. Unfortunately 
the e = 10~* contour corresponds to d < 1 pc at all 
frequencies in the LIGO band. 

Figure [41] then demonstrates that we have reached a 
second milestone not achieved in our previous all-sky 
searches [6] |7]: The multi-interferometer Hough trans- 
form search could have detected an object at the distance 
of the nearest known neutron star RX J1856.5— 3754, 



which is about 110-170 pc from Earth [Ulll^]. We could 
not have detected that particular star, since the recently 
observed 7 s rotation period [13] puts the gravitational 
wave frequency well out of the LIGO band. But the top 
of Fig. [41] shows that we could have detected a Crab- like 
pulsar (/ « 100 Hz, / « 10"^° Hz s^^) at that distance 
if gravitational radiation dominated its spin-down. 

For the ongoing S5 data run, expected to finish data 
collection in late 2007, several refinements of these meth- 
ods are under development. The StackSlide and Hough 
methods can be made more sensitive than PowerFlux by 
starting with the maximum likelihood statistic (known 
as the .^-statistic [6l HO] Il8|) rather than SFT power. 
This increases the time-baseline of the coherent step in 
a hierarchical search, though at increased computational 
cost. The lower computational cost of the Hough search 
would be an advantage in this case. Multi-interferometer 
searches also increase the sensitivity, while reducing out- 
liers (false-alarms), without having to increase greatly 
the size of the parameter space used, as illustrated by the 
Hough search in this paper. A multi-interferometer ver- 
sion of PowerFlux is under development, as well as hier- 
archical multi-interferometer searches that use the Hough 
and StackSlide method on the .F-statistic. 

Thus, PowerFlux will be the primary tool used for 
semi-coherent searches using SFTs, while the Hough and 
StackSlide methods will be used in multi-interferometer 
hierarchical searches. Strong candidates from the Pow- 
erFlux search will be fed into the latter type of search 
as well. The parameter space searches described here do 
not take into account the correlations that exist between 
points in the four or five dimensional parameter space 
(including those on the sky). A map of the mismatch be- 
tween a signal and the parameter-space templates can be 
used to generate a parameter-space metric to reduce fur- 
ther the number of points needed to conduct a search, a 
method under development for the hierarchical searches. 
Finally, the strain noise of the S5 data is lower by about 
a factor of 2, and the run will accumulate at least 1 year 
of science mode data. 



APPENDIX A: POWERFLUX POLARIZATION 
PROJECTION RELATIONS 



The PowerFlux method uses circular and four linear 
polarization "projections" to increase sensitivity to dif- 
ferent source polarizations [44|. The projections are nec- 
essarily imperfect because the interferometer itself is a 
polarimeter continually changing its orientation with re- 
spect to a source on the sky. There is "leakage" of one 
polarization into another's projection. In this appendix 
we present the formulae used by PowerFlux to define 
these imperfect projections and discuss corrections one 
can make for leakage in foUowup studies of candidates. 

As described in section [VD 3| the signal estimator used 
by PowerFlux for frequency bin k and projection polar- 
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FIG. 41: Range of the multi-interferometer Hough transform search for neutron stars spinning down solely due to gravitational 
radiation. This is a superposition of two contour plots. The sohd hnes (red in the color version) are contours of the maximum 
distance d at which a neutron star could be detected as a function of gravitational- wave frequency / and its derivative /. The 
dashed lines are contours of the corresponding ellipticity e(/, /). In concert these quantities tell us the maximum range of the 
search in terms of various populations (see text for details). 



ization angle V'' is 

where W, = is the weight for SFT i and 

-P"^'([+/x]) ^^"^ antenna pattern factor for a source with 
[-1-, x] polarization with respect to a major axis of polar- 
ization angle ij:' . 

For a source of true polarization angle "0 and plus / 
cross amplitudes Aj^ and A-^, where h'j^{t) = A+ cos(cjt4- 
$) and h'y^{t) = A^ sm{ujt + 4>), the strain amplitudes 
projected onto the + and x axes for a polarization angle 
Tp' are 

/i+ = A+ cos{ujt) cos(A-i/;) 

-Ay, sin(wi) sin(AV'), (A2) 
hx — A+ cos{ujt) sin(Ai/') 

+Ax sin(cji)cos(AV'), (A3) 

where Aip = 2{ip ~ where the SFT-dependent phase 
constant <i>o been taken to be zero, for convenience, 
and where frequency variation of the source during each 



30-minute SFT interval has been neglected. Averaging 
the detectable signal power (-F'^'(+)/i-i- + i^^'(x)^x )^ over 
one SFT interval i, one obtains approximately (neglect- 
ing antenna rotation during the half-hour interval): 

(^signal) ^\[{Fl + F^,){Al+Al) 

+{Fl-F^,){Al-Al)cos{2A^) 
+2F+Fx{Al-Al)sm{2Ai;)]. (A4) 

Note that for a linearly polarized source with polar- 
ization angle ijj ~ ij/ (so that Aip — 0) and amplitude 
A^ — /iq'", Ay. — 0, one obtains 

(^signal) = lFl{h};-^-r, (A5) 

and that for a circularly polarized source of amplitude 

A+=Ay,= h§"-, 

(^signal) = liFl+F',){h^"r, (A6) 

as expected. 
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For an average of powers from many SFT's, weighted 
according to detector noise and antenna pattern via Wi, 
the expectation value of the signal estimator depends on 
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(A7) 



where rii is the expected power from noise alone, where 
(^signal"-) is assumcd to vanish (signal uncorrelated with 
noise), and where the frequency bin index k is omitted 
for simplicity. 

For a true source with parameters ^, j4-|_, and Ax , this 
expectation value can be written: 

f3,){Al + Al) 

+(l-/32)(^'+-A2)cos(2A^) 
+2p,{Al-Al)smi2A^)], (A8) 



where the correction coefficients 
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P2 



(A9) 
(AlO) 



depend implicitly on ip' through and Fx ■ 

For a linearly polarized source with polarization angle 
ip = ip' , one obtains 



1 



and for a circularly polarized source one obtains: 
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{n{i^r) + 7;{hl^'''r{i + P2). (A12) 



These formulae permit corrections for polarization 
leakage to be applied for a hypothetical source, allow- 
ing for estimation of ip, A^, and Ax from a sampling of 
polarization projection measurements. In practice, how- 
ever, the calculation of the (3 coefficients is computation- 
ally costly in an all-sky search and is disabled by default. 
Instead, upper limits on linearly polarized sources (worst- 
case pulsar inclination) are derived from the maximum 
limit over all four linear polarization projections, as de- 
scribed in section V D 3 In foUowup investigations of 



outliers, however, these formulae permit greater discrim- 
ination of candidates, now in use for PowerFlux searches 
of the data from the ongoing S5 data run. 



SFT index i) expressing the phase in a first-order Taylor 
expansion about the midpoint time, ti/2, of the interval 
used to generate an SFT, we can write 



(j){t) = 01/2 -I- 27r/i/2(i - ti/2) , 



(Bl) 



where 4'i/2 ^^id fi/2 are the phase and frequency at time 
ti/2. Treating the values of P+ and Fx as constants equal 
to their values at time ti/2, the signal strain at discrete 
time tj is approximately. 



hj = P+A+cos(0o + 27r/i/2ij) 
+FxAxSm{(f)o + 2TTfi/2tj) , 



(B2) 



where j = gives the start time of the SFT, and (pQ is 
the approximate phase at the start of the SFT (not the 
initial phase at the start of the observation), i.e.. 



<t>o = 01/2 - 27r/i/2(r,oh/2) . 



(B3) 



Using these approximations, the Discrete Fourier Trans- 
form, given by Eq. (12 1, of hj is 
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/i/2Pcoh is usually not an 
Equation holds for < k < M/2 and 



where Ak 
integer. 

|k — fc| << M, which is true for all of the frequencies 
over which we search. 

If the discrete time samples of the data from the de- 
tector consist of a signal plus noise the expected value of 
P is approximated by 



(B5) 



where the mean value of Pq is 1 and its standard deviation 
is \/\/N due to the normalization used, and 



(rf^) = 



2 /PjsinVAK) 
+ \ Sk 7r2AK2 



Fl sin2(7rAK) 
7r2AK;2 



Pco.: 
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APPENDIX B: STACKSLIDE POWER AND 
STATISTICS 

a. Approximate Form For The StackSlide Power 

It is useful to have an analytic approximation for the 
StackSlide Power P. For a single SFT (dropping the 



is an approximate form for the square of the optimal SNR 
defined in Eq. (71) in reference [18^ averaged over SFTs 
(i.e., the angle brackets on {cP) represent an average over 
SFTs) and where for each SFT the index k is the nearest 
integer value to n. Thus, the relevant range for Ak is 

to 0.5, corresponding to a frequency mismatch of to 

1 /2 of an SFT frequency bin. 
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b. StackSlide Statistics 



III. ACKNOWLEDGMENTS 



It can be seen from Eq. (16 1 that, for Gaussian noise 
in the absence of a signal, 2NP is a variable with 2N 
degrees of freedom |TS]. Thus, the quantity 

g = 2NP 

follows the distribution: 

V{g;N)dg= ^^vw^n g^~^e~g/^dg. 



2^r(iv) 



(B7) 



(B8) 



When a signal is present, g follows a non-central distri- 
bution with 2N degrees of freedom and a non-centrality 
parameter N{(P) such that 

Vig;N{d'^))dg = 



iN{d^)) 



^-(e+N{d'))/2 



dg, (B9) 



where the form given here is based on that given in [28] , 
and In-1 is the modified Bessel function of the first kind 
and order — 1 . 



The distribution described by Eqs. (B8) and (B9) can 



be used to find the minimum optimal signal-to-noise ra- 
tio that can be detected using the StackSlide search for 
fixed false alarm and false dismissal rates, for a targeted 
search. For a 1% false alarm rate, a 10% false dis- 



missal rate, and large N Eqs. (B7) and (B9l giv e (d^ ) — 
7.3S5/^/N (See also [H]), while averaging Eq. (B6l in- 
dependently over the source sky position, inclination an- 
gle, polarization angle, and mismatch in frequency gives 
= 0.7737(4/25)(/igrcoh/5') (see also Eq. 5.35 in [10] 
). Equating these and solving for /ig, the characteris- 
tic amplitude for a targeted StackSlide search with a 1% 
false-alarm rate, 10% false-dismissal rate is: 



(ho] 



targeted 



7.7%/5/(rcohT:bJ 



1/4 



(BIO) 



where T; 



obs 



iVr^Qi, is the actual duration of the data, 
which is shorter than the total observation time, Tobs, 
because gaps exist in the data for times when the detec- 
tors were not operating in science mode. Comparing this 
expression with Eq. 5.35 in [TD] the StackSlide character- 
istic amplitude given in Eq. (BIO) is found to be about 



10% lower than a similar estimate for the standard Hough 
search. Note that in this paper an improved version of 
the Hough method is presented. Also, in this paper an 
all-sky search for the loudest StackSlide Power is car- 
ried out, covering up to 1.88 x 10® templates, and only 
the loudest StackSlide Power is returned from the search, 
corresponding to a false alarm rate of 5.32 x 10~^°. Fur- 
thermore, the upper limits are found by injecting a family 
of signals, each of which has a StackSlide Power drawn 
from a different noncentral chi-squared distribution. Us- 
ing the results from Sec. |VII[ for an all-sky StackSlide 
search the 95% confidence all-sky upper limits are found 
empirically to be approximately given by: 



(ho) 



all— sky 



23V^/(r,ohro*bs)'/'. 



(Bll) 



The authors gratefully acknowledge the support of the 
United States National Science Foundation for the con- 
struction and operation of the LIGO Laboratory and the 
Particle Physics and Astronomy Research Council of the 
United Kingdom, the Max-Planck-Society and the State 
of Niedersachsen/Germany for support of the construc- 
tion and operation of the GEO600 detector. The authors 
also gratefully acknowledge the support of the research by 
these agencies and by the Australian Research Council, 
the Natural Sciences and Engineering Research Coun- 
cil of Canada, the Council of Scientific and Industrial 
Research of India, the Department of Science and Tech- 
nology of India, the Spanish Ministerio de Educacion y 
Ciencia, The National Aeronautics and Space Adminis- 
tration, the John Simon Guggenheim Foundation, the 
Alexander von Humboldt Foundation, the Leverhulme 
Trust, the David and Lucile Packard Foundation, the Re- 
search Corporation, and the Alfred P. Sloan Foundation. 
This document has been assigned LIGO Laboratory doc- 
ument number LIGO-P060010-06-Z. 



38 



A. Abramovici et al, Science 256, 325 (1992). 

B. Barish and R. Weiss, Phys. Today 52, 44 (1999). [23 
B. Abbott et o/.(The LIGO Scientific Collaboration), [24 
Phys. Rev. D 69 102001 (2004). 

B. Abbott et al.(Thc LIGO Scientific Collaboration), M. [25 
Kramer, and A. G. Lync, Phys. Rev. Lett. 94 181103 

(2005) . 

B. Abbott et a/.(The LIGO Scientific Collaboration), M. [26 
Kramer, and A. G. Lyne, Phys Rev. D 76, 042001 (2007). [27 
B. Abbott et a/. (The LIGO Scientific Collaboration), to 
appear in Phys. Rev. D, gr-qc/0605028 (2006). 
B. Abbott ct al. (The LIGO Scientific Collaboration), 
Phys. Rev. D 72, 102004 (2005). [28 
The Einstcin@Homc project is built upon the BOINC 
(Berkeley Open Infrastructure for Network Computing) [29 
architecture described at http://boinc.berkeley.eclu/. 
Results from the distributed computing project Ein- [30 
stein@Home can be found at 

http : / / einstein.phys .uwm. edu/. [31 

B. Krishnan, A.M. Sintcs, M.A. Papa, B.F. Schutz, S. 
Frasca, and C. Palomba, Phys.Rev. D 70, 082001 (2004). [32 
M.A. Papa, B.F. Schutz, A.M. Sintes, in Gravitational [33 
waves: A challenge to theoretical astrophysics, ICTP Lec- 
ture Notes Series, Vol. Ill, edited by V. Ferrari, J.C. [34 
Miller, L. RezzoUa (Italy 2001) p. 431. 

P. Brady, T. Creighton, Phys.Rev. D 61, 082001 (2000) 

C. Cutler, I. Gholami, and B. Krishnan, Phys.Rev. D 72, 
042004 (2005). [35 
P. Brady, T. Creighton, C. Cutler and B.F. Schutz, Phys. 
Rev. D 57, 2101 (1998). [36 
G. Mendell and M. Landry, "StackSlide and Hough 
Search SNR and Statistics", LIGO technical document 
LIGO-T050003 (2005), available in [37 
http : / / admdbsrv . ligo . caltech . edu/ dec/. [38 

C. Palomba, P. Astone, S. Frasca, Class. Quant. Grav. [39 
22, S1255 (2005). 

V. Dergachev, "Description of PowerFlux Algorithms [40 
and Implementation" , LIGO technical document LIGO- 
T050186 (2005), available in [41 
http : / / admdbsrv . ligo . caltech . edu/dcc/. 
P. Jaranowski, A. Krolak, and B.F. Schutz, Phys. Rev. [42 
D 58, (1998) 063001. 

D. Sigg (for the LSC), Class. Quant. Grav. 23, S51 [43 

(2006) . 

LSC Algorithms and LALapps Applications software li- [44 
braries available at 

http: //www. Isc-group .phys .uwm. edu/daswg/. 

C. Palomba, Mon. Not. R. Astron. Soc. 359, 1150 (2005). 
R. N. Manchester, G. B. Hobbs, A. Teoh and 
M. Hobbs, Astron. J. 129, 1993 (2005). See also 



http : / /www. atnf . csiro . au/research/pulsar/psrcat/. 
S.D. Mohanty, Class. Quantum. Grav. 19, 1513 (2002). 
S.D. Mohanty, S. Mukherjee Class. Quantum. Grav. 19, 
1471 (2002). 

B. Krishnan, Bias in the estimator of the median, LIGO 
technical document T040144 (2004), available in 
http : / /admdbsrv . ligo . caltech . edu/dcc/. 
X. Siemens et al. Class. Quant. Grav. 21, S1723 (2004). 

A. Dietz et al, "Calibration of the LIGO Detectors for 
S4", LIGO technical document LIGO-T050262 (2005), 
available in 

http : //admdbsrv . ligo . caltech. edu/dcc/. 

P. Jaranowski, A. Krolak, Phys. Rev. D 61 062001 

(2000). 

The condor package is available at 
http : / /www . cs . wise . edu/condor/. 
The Matlab program is available at 
http : / /www .mathworks . com. 

P.V.C. Hough, In International Conference on High En- 
ergy Accelerators and Instrumentation, CERN (1959). 
P.V.C. Hough, U. S. Patent 3,069,654, 1962. 
J. Illingworth and J. Kittler, Computer Vision, Graphics, 
and Image Processing 44, 87 (1988). 

B. Krishnan and A. M. Sintcs, Hough search with im- 
proved sensitivity, LIGO technical document T070124 
(2007), available in 

http : / /admdbsrv . ligo . caltech . edu/dcc/. 

G. J. Feldman and R. D. Cousins, Phys.Rev. D57, 3873 

(1998). 

See EPAPS Document No. [number will be inserted by 

publisher] for numerical values of upper limits derived for 

each method in 0.25-Hz bands in the range 50-1000 Hz. 

B. J. Owen, Phys. Rev. Lett. 95, 211101 (2005). 

R. X. Xu, Astrophys. J. 596, L59 (2003). 

M. Mannarelli, K. Rajagopal and R. Sharma, arXivihep- 

ph/0702021. 

G. Ushomirsky, C. Cutler and L. Bildsten, Mon. Not. 
Roy. Astron. Soc. 319, 902 (2000). 

F. M. Walter and J. Lattimer, Astrophys. J. 576, L145 
(2002). 

M. H. van Kerkwijk and D. L. Kaplan, Astrophys. Space 
Sci. 308, 191 (2007). 

A. Tiengo and S. Mereghetti, Astrophys. J. 657, LlOl 

(2007). 

V. Dergachev and K. Riles, "PowerFlux Polarization 
Analysis", (LIGO technical document, LIGO-T050187 
(2005), available in 

http : //admdbsrv . ligo . caltech . edu/dcc/. 



39 



